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ABSTRACT 

The objective of the Eeactor Dynamics. Mcdule, ED- 1 , 
^s to obtain the kinetics equation without feedback and sclve*the 
kinetics equations numerically for .cne to six t delay €d * neutrci) groups 
for time varying reactivity • inser-ticns. The computer code • F DMCKI 
(Fundamental Mode Kinetics) wiU calculate the power- as a function of 
tim'e fcr either uraniuin or plutonium. Either fuel c^n be used with 
one to six delayed neutron groups and ona of three-types-cf. 
reactivity insertions:^ constant reactivity, sinus„c$dal, or a ramp. 
The code does not compute any parameters so the neuti;c'n generation > 
time must be provided. The user has the* option of studying the 
effects of various e time steps in solving the systenu Ihe objective of 
Module ;BD-2 is to examine the temperature feedback mechanism of a 
pressurized water reactor £P#R).and solve the one delayed neutron ■ 
ifio4el tfith temperature feedback for a* step ^insertion and a ramp 
insertion of reactivity. A* PWE Gore with -a two-path feedb.ack is 
considered. The gcre /region is the cnly one of interest in this 
module. The program ' name' is F0.MOTEM "(Fundamental Mode Kinetics with • 
Temperature Feedback) «. There are four types of reactivity inputs that 
t-fce program can accommodate* (Author)' *- 
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KINETICS MODULE '1 



THE REACTOR KINEXI6S EQUATIONS 




1.1 .Object of Modiile 

The object of this module is to: - 

• * • ' \ 

1) Obtain the Kinetics Equation without feedback 

and " 2). Solve the kinetics equations numerically for one to six delayed 

/ neutron groups forj^ime varying reactivity insertions. 

» V"* % 

The time dependence of a modern reactor is really very complicated. The 
control rod motion i*s a local perturbation soothe time dependence of the flux 
cannot/ completely be divorced from the space dependence. o The- fundamental 
mode kinetics equations do provide a rough idea of what the time behavior of 
a reactor. will be. In this module we'will develop the -kinetics equations 
and indicate How they can be solved numerically. Feedback effects will h, 
-introduced in later modules. n 

The computerize FUMOKI (Fundamental^ Mode Kinetics) will- calculate the 
power as a function of time for either uranium or pldtonium. Either fuel can 
be used with 'one to six delayed neutron groups and one of three types of 
-reactivity insertions:, . -* * 

1) a constant reactivity p< 

I J o 4 

2) ' sinusoidal p(t) * P Q sin b 2 t 

or 3) a ramp ^ m p(t) » p Q (l + b 3 t ) • 

The ramp is to simulate the rod withdraw 1 or insertion reactivity input. 

-The code^ does not confute any parameters so the neutron generation time 
inus.t be provided. Also the user has. the option of studying the effects of 
various time stepson solving the^system. The time step can'be surprisingly 
large (6.05 sec) aad stilTy^eld gob^, results in most cas^s . 1 
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1.2 The Kinetics Equations 



The time behavior of a reactor is a very important consideration in 

( 

the operation of a nuclear power plant. Also, the safety analysis o£ a 
plant depends' upon a thorough knowledge of the kinetics equations. The 
many types of reactor designs -necessitates' consideration of various 
reactivity coefficients and dynamic response characteristics. 

The neutrpnic-considerations of a reactor cannot ^be divorced from 
.associated feedback mechanisms sudr as heat transport, fluid flow, mechani- 
cal changes etc. There are many ways to delineate the dynamics problems 
of a power plant but a natural one .seems -to t?e to classify problems 
according to the , time constants involved.. There are four areas that we 

can study; •. ■» • ,* 

1) Vej^ slow transients^- fuel depletion with time constants ot 
ChB^order of a year or so. We will not consider this as a 
dynamics, problem at_all but rather in the'stati^s sections. 

2) Slow^ansients* - Xenon and slmarium effects. The time constants 
r here are of "the order of hq^irs. 

3) Normal transients - Changes in fuel and Moderator temperature, 
void changes, delayed fie utrpn considerations etc. The time 
constants here are^of /ffi^ #r*ler of a second or so and small , 
reactivity changes ar^ ^itvolved. ■ • ' * 

• . 4)- Fast transients - Control rod dropped in or withdrawn at its • 
'maximum rate etc. Reactivity inputs are of 50<? and up.. The 

time constants are of the order of 10 • sec. and serious safety 

' - r\) 
questions are raised. | 
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V 

In addition to these rather general' time reference frames, another 
important aspect of the reactor dynamics problem is whether or not the . 
core is so large that spatially dependent analysis is necessary. The 
^solution of one or more dimensional kinetics prdSlems necessitates^ the 
use of a relatively large computer. Xenon oscillations are of importance 
hefe. Also reactivity insertions are usually localized so that hot spots 
may develop. , - 

Coupled core kinetics considerations are of interest In some types 
of reactors with large cores or 'core regions which are loosely coupled. * 
Wt\en the neutron fligjit times between ^different regions of the jreactor 
are not negligible, then coupled <!ore kinetics may be a useful -tool-. 
Coupled core kinetics equations are generally differential-integral 
equations with time-lag kernels. Coupled core kinetics involves* writing' 
the kinetics eqiiations for each region of the reactor and then coupling 
the regions by neutron leakage from one region to the other. 

When a control rod in a reactor is removed, the neutron flux distri- 
bution is disturbed so that not only the fundamental mode is present but 
the higher modes are also present. However, these high^rNijarmonics die 
out very rapidly so that while the rod is in motion (in the order of seconds) 
the fundamental mode is the only important one. Therefore we will deal 
only with the fundamental mode for the space dependence and the kinetics 
equations we solve will simply indicate how the amplitude, of this £<&da- 
mental mode changes with ^tirne. For. efemple , in a spherical reactor the 
flux is • . ~ \' . % • 

• Kr.t) -V.<t>'S*§J* 



and this spatial dependence persipts through most transients, and the 
kinetics equatiorts simply yield the <|>o(t). In a word, the spatial and time 
depeniience are separable for most transients. 

¥ 

We will derive the kinetics equations using orfe group diffusiQn theory 
but the same equations are obtained 'from multi-group or transport theory. 
The thermal neutron diffusion equatio^ is 



D V 2 * - l a * + S = I|± . * (1,2.1) 



The^source is composed of three terms, i.e. 



S « S % + + S ' 

P D ext , (1.2.2) 



i 



where S = prompt neutron source, . - 

= delayed'neutron source ,> 

and S = a stfurce^bf neutrons entered externally, > as from a 
ext . 



Blutonium-Beryllium source. 
/'Not -all nei 



eutrons are emitted immediately from the fission process.' 
Th^ prompt neutrons are emitted within 10 ^ sec. of the fission process 
itself bdjt there are other neutrons called delayed neutrons. They arise 
from the beta decay of some of the fission fragments. As an example, 
Kr-93 is a fission fragment which has a half .life of 1.22 sec. This beta- 
decays to Rb-93, but it is formed in such a hifehly excited state that a 



neutron is emitted following the beta decay and the reaction is thus, 



■ ■ xi,.c. »• 'V— 35 — ► ■ 

"10. sec. o 
93 

There are some 30 of these neutron precursors^ Kr) which produce neutrons 
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fhrough beta decay. ^ 

r • > 
* Since the half-lives of many of , these precursors are ver^clo^e to 

each other, it is necessary only to consider six delayed neutron precursor 

groyps which are averages of the 30 or so precursors, appropriately 

weighted. -The delayed neutron precursors we talk about are thus fictitious 

lit that they do not actually exist, but are averages of the afctudl* 

precursors. Table 1. 5. 1 (section 5). indicates the groups of delayed 

\ neutrons obtained from .these six precufsprs for neutron induced fission 

4ji uranium and plutonium. 



• 235 
The fraction 1 of neutrons which are delayed is called B. For U, 



B^'o.0064 and the total delayed neutron fraction is .the sum of the 
individual pnes , i.e., r • 

i-1 , • • 
Returning now to the diffusion equation, we can, set 

S = vL *(r, t) . • (1-2.4) 

" -* ? f • 

/♦ n (1.2.5) 
S D | H A C (r, t) . . 

" a i-1 » 

where C i4> (r, t) is the precursor concentration of group "i", i.e., the 
numbers pre cursors /cm 3 existing at 'point r at tWt.~V Each time a 
precursor C dec^s, it is assumed a 'delayed neutron of group "i" results. 
So our diffusidh equation becomes (Equation 1.2.1) • 




a «v * i*l" 



(1.2.6) 



r 
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Equation(1.2.6)has four independent variables, x, y, z and t, and' 
seven dependent variables C^, C 2> . . . , v C fi . In order to solve such a 
system of equations, we will have^to write six more equations, one for- 
each \ of the precursor concentrations. The rate at which the precursor 
, concentration changes is 

3 C. (?, t) - 6 ± vE f <K?, t) - A i Q ± (r, t). - (1.2.7) 



3t — : 1 L 



rate of formation rate of decay . 

i - 1,,2, ...,6. 

This system of equations is relatively difficult to solve so we shall 

work analytically with Equation (1.2.6) and (1.2.7) for awhile.* 

^ Assume that we want to expand the spatial dependence of the flux 

$(r, t) and the precursor concentrations (r, t) in terms of a set of 

» 

v. 

eigenfunctions of the Helnholtz equation 

V 2 Y (r) + B 2 Y (r) *'0. - (1.2.8) 

n n n 

The reason for doing this is that the solution *°r thfe steady state satis- 
fie s/ this same equation with coefficients which are time independent. In 
the time dbp^ndent problems the coefficients are functions of time.. The 

"effcfen functions Y\*(r) are Cos B x for a slab reactor, Bessel's 

- ) * n n 

functions for a cylindrical reactor, etc., -In' any event, we set 

CO * ^ > 

* (r, t) - E 6 (t) Y (?) • .(1.2.9) 

n-0 . . 



J 



k 



and 

v 



00 



C' (r/t) - E C. (t).Y (r). (1.2.10) 
n=0 



We note" that <t> n (t) will yield the^amplit\id#of the ti~ harmonic -and 
Y n (r) its spatial distribution. Substituting these ejxpressions into 
Equations (1.2.6) and (1.2.7) and using Equation (1.2.8) we obtain 

r 

E [>db 2 *_(t) - z a i n (t) + (l-e) vz f ,J n (t) |Y n (?) 

„n=0 L • — 1 _ 



i-l „, 0 , 1 " - n-0 ". it 



and 



£ Y n <*> '^aT^ C 6 i VZ f *n (t) "' X i. a in <£| *»•<*>■ 

•n=0 n=0 



^Now using the fact that the Y (r) form an orthogonal set, the preceding 
equation is simplified by multiplying it by Y m (r) and integrating over 
the reactor volume (taking the inner product) to yield • .> ' 

mm * 

1 d<t> n " (t) = -(DB 2 + I -M (t) + (1-3)' vl f * (t) + 5 ■ 
— /- — nan r n . 

v dt . 

6 s . ' ' . 
E X. C in (t),« . , (1.2.11)' 
i-l . 



and .9 • ■ ~j 

dC, (t) - ' 

\ = UU n (t) - . X. C in (t). • (1.2.12) 

y^x -jp i £ n . i , - 

V 

Nt>w it is necessary to solve these *wo equations for the expansion 

coefficients * (t) and C (t) . .This is a very difficult task so we 
t in . \ 



i2 •. 



make some definitions and some approximations. 
* First we define the multiplication k as 



. * 4 



v Zjl ' J - 
f a » 



k be 1— S— • (1.2.13) 

n 1 + B 2 L 2 
n 

: - 
and the thejmal neutron lifetime ' ^ 



it 



n , v E- (1 + B 2 .L 2 ) • ' (1.2.14) 

an 

The next thing we do is assume '^hat there are no delayed neutrons (there 
are, but assume for the moment we can neglect them). Then Equation (1.2.11) 
becomes (if 8 * 0) , 

d K (t) = -v I (1 + L 2 B ) 4 (t) + v VE 4 (t). , ; T 

- ■■ " - a n n in 

dt ^ 



Now using Equation (1,2.14) we have 



) 



- — L 2* an, 

dt 4 au a 

n 

« 

2 .2 

4 (t) 1 + B n L • 

* a " 1 + B 2 L 2 n 

11 ■ n 



n « ^ 

The solytion of Equation (1.2.15) is obviously 

• V 1 

• - , x "T~ C (1.2.16) 

4 (t) « 4 (o) e n 

-n n 
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^For^a typical light, water reactor, l Q - I - 10~ sec. and k Q - k • ,1.0020 
From Equatiqn Uy2.14) we note that * v , 6 * 



2 T 2 -2 



r 



f 




so that the higher harmonics die out rapidly. We will keep only the 
• • # * < 

fundamental mode, not only 'in the solution of the no delayed neutron case 
' but also for the solution of Equations (1.2*11) and (1.2.12).- 
w Returning now to Equ^ions (1.2,11) and (1.2.12) and omitting the 
•subscript n since we are 'concerned only with* the fundamental mode, we s«t 



, ' <J> (t) » v n (t) 8 



-3 

where n (t) is the neutron ^density (cm ), apd after a bit of algebra 

j 

come up with' * ^ , * 4 4 

dn(t) =? (1-B) V" v E n (t) - v E " (1 + B 2 L 2 ) n (tV + 21 X C.(t) § 
dt * ' a a . 3 i-1 1 



5* «** • 

and • 

d C (t) Z 

• .V * "it = 8 i vZ a v r n(t) " A i C i (t) ' 

^ v a 



'ffc*. Using the definitions given in Equations (1.2.13) and (1.2.14) the above 
equations become 




i=i 
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k'is the effective multiplication, and the effective neutron lifetime is 
depbted as t. If we drop thV tilde on the C^s', then we have > ' 

° r - Mtl' J iV- 6) kn'l B (t) ; C. (t) J^.2.18) 



and 



dC i (t) ■ = 6. 7 'n(t) -X. C. (t). (1.2.JL9) 



These are* called the fundamental mode -reactor kinetics equations, in 
that spatial and time dependence are assumed ^separable* 

Sometimes it's convenient to cast these equations into a slightly 
different form. The thermal neutron lifetime is 

X - 

i «= -& V-9 4 * - • (1.2.20) 

V • 1 + B 2 L 2 ' 

i.e., the mean time a neutron spends in'the system from its birth as a 
thermal neutron until it's absorbed or w leaks out of the reactor. The 
neutron gereratiqji timfe is defined as- 



\ 2*2 

V . ' . 1(1 + L ) • t 

A >.£ - ■ 1 TT-^T . - ; = TTTF" ' (1.2.21) 



' k - vl .(1+B 2 L?) ' VE f , VVE f 

a 



- v - . • 



The "generation time is the|pin time that it takes one neutron to generate 

one more -prompt- neutron or\one precursor. The neutron lifetime is thus 

the reciprocal of the .destruction rate of neutrons and the generation 

» » 

time is the reciprocal of the production rate of neutrons. With this 
definition and the definition of reactivity p n 



pH'.f , • (1-2.22) 



t 

the kinetics equations (1.2.18 and 1.2.19) become ; 

dh_(t) = p(t)- - B , ( j + y x c (t) (1.2.23) 
dt * fci 1 1 

d c i (t> . *-± nCt) - • X, C. Ct) (1.2.24) 

9 * -.A 11 , 

, _ ( i = 1,2,. ...6 ' 1 

Equations' (l.*2.*23) and (1.2.24) are^used to describe the time behavior 
of a reactor. . This form <Jf the equations results if multigroup diffusion 
theory or transport theory is used to derive the. equations but the 

definitions" of £, A, p, k are modified. WeH/ill assume they are- input 

* * - 

parameters to the^ystem of equations. _ • 



7* 
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Problem 1.2.1 



Estimate the neutron lifetime and neutron generation time in an 

k * y 

infinite stack of graphite and U-235 if a «- 5 x ICT b for the graphite 

a . 

and there are 500 atoms o^carbon to each atqm of U-235^ There is no 
U~238-;present. o « 680 b for the U-235 and, take v to be 2200 m/sec. 

/ • v-. ' a ./ " v 

Problem 1.2.2 . ' ' 

Estimate the pferiod *of^ a reactor if there are no delayed neutrons 

f 4 v 
-4 . > 

and if the k = 1.0010 and i =* 10 sec. The reactor period is the time . 

it takes for the flux or neutron density to increase by a factor of e.* 



Problem 1,2.3 

A X 



i 



Let Y_ L » C ± (t) -g — and = g/A. Show that the kinetics 
equations are then 



* ft 



dn 
dt 



[n(t) (p»-l); + T a i Y ± ] 



i-1 



and 



j~ - .*i [»CO - x 1 (t) t ] 



Where p • ■= arid a 4 = 
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r 

* 0 



t 

1.3 Analytical Solutions of the Reactor Kine tics Equations 

o . « v * 

The kinetics equations are relatively difficult to soJveLboth 
'analytically and numerically due to the large difference in the time' constants 
in the equations'"^ well as due to the fact that ttyjre are seven coupled 
ordinary differential equations if there are six delayed neutron groups. 
Equations (1.2.23) and (1.2.24) are repeated here as 



i=l 



(1.3.1) 



and 



dC.(t) 
~dt 



^i n(t) 



i - 1, 2 6. 



(1.3.2) 



If we assume' that we can take an appropriate average for'the delayed 

neutron decay constant, then it is possib-le to collapse these seven s 

* > ' 
equations ,into two equations. For one effective, group of delayed neutrons, 

t • ' ♦ 7 

EquationsJ>1.3. 1) and (1.3.2) reduce to* • 



' • / Mt) ," f(t) - 8 n(t) + \ C (t) 

,( dt A 



(1.3.3) 



and 



dC(t) 
dt 



h(t) - XC(t) 



(1.3.4) 



The ""decay constant. A for this average delayed neutron precursor C(t) is 



1 1 f ii 
X " 6 £i X i 



(1.3.5) . 
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This average is not unique and at times it may be .advantageous to use 
another expression for the average/ Tf very small reactivities are added 
to the'reactorl then instead of .averaging over all six groups, perhaps 
only the longest three would be used. 

If only one group is desiijed Equation (1.3.5) is used to determine the 
effective A and $ is simply the sum of the 8^ For two groups, we split the 
g's into. two groups of three each and 



1 i=l 1 X i B l 1=1 \ 

- m 



(1.3.5a) 



■J 



' and 



6 ' • V 6 3 ' 

= T 8 - - T • (1.3.5b) 

& X 2 ^2 l " 4 X i . • 



For three groups, the 3 ! s are split into three groups of two components each 

and for four groups, .the shortest two groups are averaged together as well 

as the next shortest, two and the^22 and 55 sec!, group? are treated separately. 

We will solve Equations V (1.3.3) and (1.3.4) subject to a constant ^ 
reactivity insertion P q at time zero*. In ordfer to do this we recast otjr 
equations into a matrix form because of the similarity to the numerical 
i technique (Hansen's method) us.ed in the solution of the equations. In matrix 
form, Equations (1.3.3) and <1.3,4) become ' « 

di(t) 

l -j^- - A± (t) (1 . 3f6) 



with 



±(t) 



n (t) 



C (t) 



V 



and A 



P o -fi 



£ 
A 



Notice that the A matrix is independent of time if the reactivity is a 



constant. Equation (1.3.6) can formally be integrated to yield 

± (t) = e = . i (0) . 



(1.3.7) 



The initial condition vector ^ (0) is 



i (0) 



"n(0)" 


- n(0)" 


" 1 


C(0) 




B 



(1.3.8) 



where we have used Equation (1.3.4) and set the derivative — = 0 for t < 
WhenOreactiv^ty is added at time zero, the delayed neutron precursor 
concentration C(t) does not change until after some time. It is assumed 
that the, reactor has been operating for a jLong time at a reactor power 
consistent wijjjff a neutron density nCO)'. 

The formal solution given by Equation '(1.3. 7) is not much gpod if the 
explicit functional dependence of n(t) and C(t) cannot be obtained. The 
reason po explicit functional dependence is achieved £rom Equation' (1, 3. 7) 
is that matrix A 'is not diagonal. This means we really haven' t separated 
the equations from .each other. In order to demonstrate how. the solution 
can be obtained, we imitate the teckoiq ucTused for the solution of an 




v- ■ / 

n— order differential equation, i.e. §et 



i(t) ^ v , 



/ 

(1,3.9) 



where v is a constant .vector .and 0) is a scalar independent of time. 
Substituting this into Equation (1.3.6) we have ' " f 



1 0)t \ ■ o)t . 

0) e y . =■ A e v 



or 



A v = w v; 



(1.3.10) 



Now it is apparent that Equation (1. 3.10) is an eigenvalue equation so 



I A" - tulft '= 0 



9 



(1.3.11) 



jnust be satisfied to determine the eigenvalue to." This equation is Called 
the characteristic equation and in reactor dynamics it is called 0 the H in- 

$ . . . o j * O 

hour equation" because it .relates the reactivity 0 inserted into A with 
the oj which is the reciprocal ql * the '-reactor period. ? 

The inhour equation (Equation 1.3.11) for the one group delayed neutron 
model is , - 



det 



A 



- 0) 



— OC1..3..12) 



or 



*Aw+(0-p< + Aui)-Xp ^0.; 



V 



(1.3.13) 



Equation (1.3.13) can also be written as. 



A -a) +' 



(1.3.14). 



If we had* k£pt all six groups of delayed neutrons, Equation (1.3.11) 
rwould be a 7 x 7 determinant and Equation (1.3.14) would be enlarged to 

\ 



' ' 6 Sm 
i»l i 




(1.3.15) 



The v of Equation (1,3.10) is'thfe eigenvector associated with the 
eigenvalue u). .Now if we perform a .similarity transformation B on A, then 
both A and B A B tiave the' same, characteristic equation. Also the fact 
the fact -that the trace of & is the^sum of the eigenvalues of the matrix A 
is a useful check on the actual' calculation of the eigenvalues of A. Frdm 
this lasfcifact it is apparent that the sum of the two roofs andP» 2 (the 
two eigenvalues) of Equation '(1.'3. 13) is given by the relation 



6 -0 
o 



(1.3.16) 



The eigenvectors' V, and v_ 2 .associated with u ]L and «) 2 respectively 
are obtained by taking the cpfactors of the element in any row of . 
Equation (1.3. 12). To see how this -works , we set * 



'11 



'21 



'.and v^ 



Nit, 

V 12 



'22 



and have. v^, the cofactor of .the element a^ ' u)^, as 



'11 



'and the cofactor of the element X of the first "row as 

v n =. - B/A . . 



Similarly, for the vector v** we have for the .second row of |A - Ico| , 



v l2 = - X , 



V 



P 0 - B 



Notic 




22 _ . A. 



ce that po get the y^, we us^the^i gen value co^ and for v 2> we use 



the .eigenvalue u> 2 in Equation (1.3.12). Therefore, the eigenvectors of 



A 



are (to within a constant) 

X + (0, 



*1 



A 



r* 



and v n 



,A 



+ to. 



These vectors v and v„ are linearly independent since the eigenvalues 
are distinct. * x> 



Now since we have 



for i = 1, 2 , 



it is apparent thJt the similarity transformation 

- .V, 

*A B ■ B D 



(1.3.17) 



hotds r where D is a diagonal 'matrix having the eigenvalues as its elements, 



23. 



■A 



i.e. 



, -- 19 - 



Wo 



and 



These two matrices are now eoitipletely determined JEor the one Relayed 
neutron model. v 



If we make a transformation 

' * 4 <t) - B Z(t) , 



(1.3.18) 



then 



><g(t) 
dt 



= B -1 A B Z(t) 



D • Z(t) . 



and its solution is obviously 



t(t) 



o> 2 t 



Dt 



Z(o) - e" . Z(o) ; (1.3.19) 



Now using- Equations (1.3.18) and '( 1.3.19) we have 



Dt 



-1 



<l(t) - J e - | i (o) . 



(1.3.20) 



( 
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This ds the solution of the kinetics equations.' Notice now that the exponential 
Dt 

e~ is a diagonal matrix- so the separation of the equations is effected. 
* ' We now write out the details of this procedure. The trans formatioQ 
matrix B is * ' 



X + a). 



B/A 



0 , 



and the determinant of B is 



det B = + — — + Xu 2 



'(1.3.21) 



The inverse % of B is 



det fl 



6-P 0 . 



+ . "2 



X + » 



The solution of the one delayed neutron group^ equations is obtained 
from Equations (1/3.20) and (1.3.8) : 



±(t) 



det B 



X + w 



- 21 - 



\ 



< 



. T"" 2 

\ . 



io 2 t 



• 1 -i 

^ -X- 



X + a), 



AX 



_ n(t) 
Det B 



p o)(X + (oj e 
C« 2 " r 1 



o) 1 R a) 2 t 

+ — e 

« A 



B-p 



0 x 1 ^ 



+ » 2 > xT c 



(1.3.23) 



Even though Equations (1.3.23) represent 1 the exact solution bf the problem, 
we are generally not concerned *&ith c(t) so we look only at the n(^t) 
equation. Also from Equation - (1. 3. 13) we have 9 



J l,2 



-(^-p + AX) ±\/(e-p + AX) 2 < + 4 XA p 



2 A 



^1.3.24> - 
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A 



Note too that 



w l W 2 = " T" * (1.3.25) 



and in agreement with Equation (1,3.16), 



X 



V G-P + AA / 

^* O^ + u^- — ^ , ; ? (1.3.26) 



Now 'if we approximate a>. (the' positive sign) and (the negative sign) 
, in Equation (1.3.24) by assuming that 

♦ 

(B - P Q + A A) 2 > > 4 A Pq A , 

i 

which is true for most reactors, we have the radical of Equation (1.3.24) 
equal to 



or 



• \ 



4A p A 

(3 - P Q + A A) ^1 ^. Pq ° + aa) 2 



|* + 2 (B-p -fr AA) z + - '."J ^P 0 + AA > '• 



When this is used we obtain 



, « • Xp B-P n 

v> • . w i 4 $4r 311(1 w 2 V 1 - ( ^- 3 - 27) 



Note 'that if p q is y positive, then.w^ is also positive. 



The determinant of § is involved in the solution of n(t) so with - 

the above* approximations as well 'as J>y using Equation (1.3.25) we have ^ 
* T ~ — — ~ ^ 



— 



* , - 23 - 




*<>„ 

o 



B _ M.. (1. 3.28) 



N o Equation U.3.23) becomes, usii)g Equation (1.3.28), 



n(t) = n(o) 



' (1.3^9) 



This approximate solution is\ery useful in obtaining checks on numerical. 

N t 
solutions. The character of thi solution is also, exponential. This will- 

* / 

play a role in our numerical/technique . 

Problem 1.3.1 , 

. — — — ^ 

' Assume that p(t) « 0 in Equations (1.3.3) and (1.3.4) and obtain a solu 



tion of the' kinetics equations. 
Problem 1.3.2 ' *" 

Prove the theorems 
a) If" 1 A B - u l| = | A - a)l| 

f 

and 

2 

- b) trace A ^ Z w, 

= i«l 1 • 



Problem ,1:3. 3 

Why does' Equation (1.3.17) follow from the eigenvector equation 

. A V e» 0) V?* 

; 

r 

Problem 1.3.4 

;| Using the* same approximations as used in developing Equation (1.3. 29) 

establish a relation, for c(t). , 

fERIC , y 7— -V 
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1 . 4 Numerical Solu^oti of the Kinetics Equations 

c 

The reactor kinetics equations are difficult to solve mu&erically 
by ''standard 11 Runge-Kutta of predictor-corrector methods. The basic 
reason becomes apparent by looking at the one-delayed neutron group 

* • y ' 

equations which we ^repeat as 



and 



f n(t> - X c(t) .• ' (1.4.2) 



The very short* time response A of the prompt neutrons -is .of the order of 
10 ^ sec. whereas the delayed neutron time response is -j~ or about 10 sec, 
a factor of 10 greater. The implication of these facts. is that in order to 
obtain the prompt response, very small time steps, of the order of 10 sec, 
are required. But then before the delayed neutron term can come into play, 
many time steps are required. Also, to examine* the response out td even one 

second, 10,000 steps of calculation would be required.. f v „ 

■ „ f * 

There are several methods that are used to ameliorate .this difficult 
problem: t 0 

1) Using a Laplace transform technique. # • 

2) Transferring the differential equations to integral equations. 

\> " 

3) . Using the eigenvalue method. . 

* * ' v 

We choose the last technique and refer ta it as Hansen 1 s method (4) after 

its ojrij|lnator. The method works for varying reactivity and can be extended 
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to systems with feedback. 

' ' ' The basic idea of Hansen's method il relatively simple. We again 
write- Equations (1.4.1) and (1.4.2) as / matrix and set 



dMt) 
dt 



A ± 



(1.4.3) 



where A and ± are [defined as in Section b of this' module. We will also 
only perform *he operations for our one flayed neutron group model. 

4 * 

Now set. 

- . a = L + D + U (1.4.4) 



where 



L -= 



r 9 
0 



rp-g 

A 



and 



T 0 ' 



0 ' 



Of course, for all six delayed neut recoups these matrices would still 
have the same meaning. Equation (1.4.4) can be written Using'these 
definitions as * 



dt 1 



-'Ji(t) - (L + U) i-(t) 



(1.4.5) 



Equation (1.4.5) does not appear any simpler to solve than Equation 
(1.4.3) and in fact it isn't. The reason for splitting it up in this fashion 
ds to develop an iteration procedure. We assume we begin this Calculation 
from a time t^ and advance to a time t^. We set , / 

•".'*(•• 

h = t, - t. . • (1.4.6) 

. -i u • 

SincV^is.a diagonal matrix, an integrating factor for Equation 

(1.4.5) is e if the reactivity doesn't change mtjich during the time interval 
h. Therefore, Equation (1.4.5) becomes 

-Dt djji(t) -Dt -Dt 

e = — — - e " D ±(t) = e 3 (L + U) jfc(t) , ' , 



or 



, -Dt 1 -Dt 



Integrating rtow from 0 to h we have 



-Dt h ' . -Dt 

e 83 Ht) I » r e 83 (l + vy±(t) dt 



or 



^ Dh H D(h-0) * 

♦(t + fc) « e S i(t ) + 'r e" - (L + U) ±{t + 0)d, (1.4.7) 

O 0$ O to « o 



where 



and obv±otisly in this interval 

dO » dt. 



Thisis^an integral equation since the" function we're looking for 
is part of the integrand. In order to provide a reasonable approximation 
to jKt Q + O) , we recall that the analytical solution is exponential so we 



assume 



± (r Q + 0) = e^ 0 i (t Q ) , U.4.8) 

f 

4 

and u is the largest eigenvalue of the matrix A. This means that We will^ 

O 4 „ 

have to solve the equation (notice this is" simply the inhour equation) % : 



A 




• |A - Hi | = o (1.4.9) 

in each time 'intervaINf or which the reactivity has changed>'since the 

/ N . 

reactivity 'will generally be a' time' dependent quantity. 

Inserting Equation (1.4.8) into (1.4.7) we have t 

Dh . g(h-Q) "8/ 

. 4<t + h)*e" l(t, 0 ) +/ o e (| + y> e i (t o )dG 

Dh r . D(h-O) u 10 "1 

= e = i(t o ) +|/ h o e- e °" dGj (L+U) ±(.t Q ) 

Dh |>. .fe (O)„I-D)0 

= e , i(t ) + 



fDh .- .6 (w I-D)© I , 

|e" • / h Q e 0= dGj i<t o ) . 

D , ru hi Dhl 

e S ^(t^+^I-D)" 1 e ° " - e J(L+U) i(t Q ) . (1.4.10) 



If we write 



¥ (t Q ) = *, 



(1.4.11) 



and 



t (t„ + h) = # 



Then Equation (1.4.10) becomes 



i 

/ 



+1 



Dh 



f id hi Dh"| 

• 0 i;S) L e ° e= J ( = + = )Jl b ; 



= G 



(1.4.11$/ 
(1.4.12) 



This G matrix obviously represents 



Dh 

G = e = + (to 



. r co' fh bh*i 



(1.4.23) 



, and it can be written as^ 



. G 



io h ,, , 

e ° - e" Xh I 

to + A, A 

o 1 • 



to h + h 



o if ' 



-Xh 



(1.4.14) 



Por completeness, if there are N delayed neutron groups, we include 



the following relation for G, 
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# h 
o 

e - e 



A . 



p-P 

o A 
e „ - e 



w h -X t h 
o 1 
e - e 

a) + A- 

o • 1 



(erg.) 

1 A ; 



e 



a) h 
o 

e - e 



-V 1 



Buz 



-v 



(1.4.15) 



If there are six delayed neutron groups , it 'is apparent that G will be 

a 7 x 7 'matrix. Also this iteration technique as given by Equation (1.4.12) 

is unconditionally stable and yields the asymptotically correct eigen- 

• * » 

solution 

* r 

The numerical procedure for the solution of the kinetics equations 
is thus: 

1. , Determine the nunber of delayed neutron groups desired and 
* read in the pertinent parameters such as A^ ^ etc. Also 

choose a time step "h". 
'* 2. Construct the vector 1(0)-. This will usually be. ^ 



£(0) = n(0) 



3. Determine the largest eigenvalue of the equation 



4 



.A * wl - 0. . 



This is a rather dif f icult~ s-tep since ifr means solving an 
algebraic equation of perhaps degree 7 to -deterjpine its largest 

root ca . J » 1 - " 

o — 

4. Construct the G matrix using Equation (1.4.15). * 

5. Determine the vector j/, where 

ijHtf) « ^ G jp(0) . 

6. Repeat the above steps starting with step 3. 
The technique is ^^R^ot involved and can yield very accurate* results*. 

The determination of the root of an algebraic/equation needs some 

* f •. * 

discussion. We qhoose the Newton-Raphson technique to solve the* equation 

N 

, =H a oi n = 0 foe, 1 KN £ 7. ' (1.4.16) 

h=0 . 

\ . • • * . ' . 
Let the N roots of Eqjuation (1.4.17) be labelled ai o> .u^ where 

^the roots are ordered such that u> q > > u> 2 < . . > uy For our problem, 
ell the roots are real and'there will|^nly be one positive .root depending 
on whether p is positive at the time step of interest. We are interested 
only in u> q . Also we assume tha^ the_interval t of interest in the roots is 
limited by - "* 



o — 1 A 1 



and 



IS 4 



with the decay cons tmt corresponding to the shortest lived delayed 
neutron group. 'The reactivity can ran^fi from negative 3. to positve 3 in 



most practical situations. Note that if 

p : » 0 then to Q =■ 0 , 

* V - 

and if P - B, then ^ - - while if p £-6, then u> 0 - -X r »«* the last _. 
Situation," regardless of the amount ofnegativ* reactivity Introduced ' 
into the reactor,' the reactor cannot shut W f*t#r than, a period of 

T a o) " 1 *» ~— " = 80 sec. 

The Newton-Raphson method is relatively simple t«*use. Assume that , 
we can expand f(») in a Taylor series about where ^ is the root of- 

interest./ Then ' 

. . . df{u~^ . * 2 d2f (% > 

f( % ) = f(. 0i ) + .h + h_ + . . . (1.4.17) 



4 , 

« is a first "guess" 'at the solution which we assume to be p(t)/A. 

„ to . m a 17 \ affpr the* firs£ two terms on the right, 
If we' truncate Equation (l.^.l/A arter cne l±l^ 



The 

o 



we have 



f ( U / - f <• > + «f ( ,V -<> 



or 



i_ (1.4.18) 



dfU ) 
du 1 



The next approximation to w o is then 




u ■ to + h 
02 oi - 



(1.4.19) 



and then Equation, (1.4.18) is used again. If this iteration procedure is 
used a sufficient number of times, the?!? the roots of Equation (1.4.16) can 
be obtained. 

Problem 1.4.1 * 

* Use the Newton-Raphson method to* solve the equation 

< 

' * : x 2 + 3x -• < C*^ 0. • • .» 

I 

Check by the quadratic formula. . , ✓ 

Problem, 1.4.2 ' • 

X • • Extend the Newton-Raphson method to systems of equations. In particular 

f • • ^ 



*v - 



if 



and 



f 1 (* v x 2 ) = 0 



f 2 (x lf x 2 ) - 0 



then show that if h^ and h 2 are the increments to the 'assumed roots, 



f 2 3jc " f l 3x 0 

v h. : 1- — L 

1 » det.J 

~ at 

and - 

• ' • ■ . 3^ 9£i 

* -* f l 3x 2 " f 2 &jc 
h 2 * ». ^ det J 

*riiere J is the Jafcobian matrix, i.e. * 4 ' „ 
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J = 



3x, 



3x, 



5 

3x„ 



and det £ is the determinant 'of the Jacobian matrix. 

Using this formulation ""find the solution to the equations 



f 1 (x 1 , x 2 ) = sin x x x 2 - x x + x 2 



2 2 

f (x v x 2 ) = 2x x - x 2 >+ xj x 2 = 0. 



1 



\ 
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1.5 The Computer, Program and the^ Kinetics Equations 

As^has already been indicated, the kinetics/ equations offer a challenge 
for their successful computer solution. The purpose of this section i3 to 
. indicate what models the program will solve. The program will solve the 
kinetics equations by means of Hansen's method for the following reactivity 
inputs : - 

1) A constant reactivity, 

P(t) - p o . • (1.5..1) 

2) Sinsusoidal variation of reactivity, 

a 

p(t) » sin b.t, (1.5.2) 

where p Q .and mu a.t be input. * j ^ 

3) Linear reactivity insertion to approximate the insertion 
> or withdrawal of a control rod: 

p(tr« p o ju + b 3 o. \ * ° (r.5.*3) 

These inputs. are-thtr~most common and can be used io simulate many 
situations. In each of the above situations p Q must be* read into the 
computer as well as the or rate of insertion in sec \ The reactivity 
units are in dollars. A dollar Of t reactivity is the amount inserted or 
withdrawn which equals the delayed neutron fraction. For example, a reactivity 

m i 

of 50<? .for U-235 vhere 0 = 0.0065 is 0.00325 whereas for Pu-239, it is 
0.00135 since^p » 0.0027. ' * 

The program also has an option of either doing the calculations for 

v - ^ • 
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235 or Pu-239. The delayed neutron decay' ccttStants are taken from Keepins' 

.* * 
data and- are^given in Table 1.5.1. 

Table 1.5.1 Delayed-Neutron Half-Live), Decay Constants and 
Yields from Thermal; F&sicm of U-235 and Pu-239. 



URANIUM-235 



* .» PLUTONIUM-239 



\. ' • 

9 


Group 
Index 


Half-Life 
(sec)- 


Decay 
Cons tan J 
, X (sec ) 


Relative 
Abundance 


Half-Life 
(sec) 


Decay 
Cons tan J 
X (sec ) 


Relative 
Abundance 
Bi/B 

\—l 




1 


' 55.72 


' 0.0124 


✓ 

0.033 


54.28 


0.0128 


0.035 

■) • 




2 


22.72 


0.0305 


0.219 


•23.04 


0.0301 


0.298 


* 


« 

3 


. 6.22 


0.3,11 


o:i96 


5.^0 


0.124 r 


0.211 


» 


4 ' 


2.30 


0.,301 


0,395 


2 A3 


0.325 


0.326 




5 * 
6 


0.6,10, 
0.230 


- 3.01 
'8=0.0065 


0.115 
0.042 


0.618 

\ 0 . 25<7 

,1 


1.12 
2.69 

8-0.0027 


0.044 



The' computer program for this module is good raly for low power reactors 
sWe feedback effects are .not taken into account. It does give some idea 
of the behavior of the' pqwer of th* reactor for various reactivities as-well 
as for the two different fuels. 



1.6 Input-Output Data for Code FUMOKI 

The- input data required for the program are presented below: ^ 

Card 1. » IFUKL indicates the type of fuel the reactor has. It can either 
235 239 

• --—be, U or • Pu. If IFUEL = 5 - its Uranium -235 / 

I FUEL = 9 - the fuel is Plutonium 239. 
The format is 12. * 
Card 2. NN - The number of groups of delayed neutrons desired. NN can 
any integer from 1 to' 6. 

The format is again 12. 

* k A , » * 
Card 3. NRO - Type of reactivity inserted into the reactor. 

' . The format is 12. \ 

If NRO = p(t) = p o , ' * 

NRO> - 2, p(t) - p n sin b„t, 

. ' , NRO • = 3i p(t) « P Q (1 + b 3 t). 

Card 3 must have the reactivity p^ , in dollars, read in with an 

F 12.9 format. The b 0 and b~ are read in with units of sec 

„ V The b's are read with an. F 12.9 format on this card. . 

Card if. XL - Neutron generation time. 

XL must be in seconds and is 'in an F 12. 9t format. 

' - • * v 

Card 5. TIME: - the total time period -for which the solution is desired 

■9 

(seconds) 

H The time step desired, again in seconds. This cannot be 
taken arbitrarily large. Generally pne should choose an 
, H of 'between ^mb^t(C^^Q^t. . Both of these are F 10.-5 ;< 
* formats. . 
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Card 6. IG - Output indicator. This indicates whether the user would 

like to print out the G matrix, the 'flux and tJjne step or ^ 
not. If IG » 1, then the flux, "G" matrix arid time 
step are all printed. If IG = 2, G is not printed out. 
This is an 12 format. 

Card 7. XN - The initial neutron density. 

XN can be anything you choose hut it is in an F 10.5 format. 

An example -of an input data set is given below. This is an example 
where we use Uranium fuel, with one group of delayed neutrons for a constant 
reactivity input of about 34<£. 

- — ~ : ' _"_ _ , " _J_ 

"zzzzr ^ , i— — h_jizh-.- ^- 

— 1 — ~"" ■ . "■ \, • # 

~ — - ▼ - y \ , 

^ • * ' - 6 ' 

' • ••X--- - . ' - " • " . 1 ' 



2 2 2 2\,2 2 2 2 2 2 2 ^ 2 2 2 2 2 7 2 2 2 2 2 2 7 2 2 2 2^ 2 '22 2 ? 2 2 1 i 2 4 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 Z 2 2 2.2 2 2 2 2 2-J 2 2 2 2 2 2-2 2 2 ^ 
33 3 3*33 3335 3 3 33 3 3 333 33 3333 3 3.3 3 3 333 3 3 33 nlgi 3 3 3 3 3*3*333333333 3 3333 3 3 33 31333 3 3 3.3 3J* 



44444444 4* 4444444444444 4. 4 4444 4 4 444444 M 44 

54.5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 55 5 5 5 5 >56 5 5 S ; 
6,6 6 6 6 S 6 6 6 6 6 G G G 6 6 6 6 6 6 6 6 G 5 6 6 6 & 6 G 6 5 G G 
M 7 7 1 1 7 7 1 1 1 11 1 1 1 1 1 1 1 7 7 1 1 7 7 7 7 7 7 7 7 7 7 7* 
8 8 8 8 ? 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8 8*8 « 8 8 $ t f I * • • 




4444444 4 444 4 4 4 4 444 4 4*4 444 4 4 4 4 444 4 

5 5 5 5 S5 5 5 5 5 5 1» 5 5 5555^5555555555 5.S 59 
G 6 6 6 8 G 6 G G 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 b 6 6'6 6 6 6 8 6 



M*l 7 7 ? 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 T 7 7 7 ? 7 7 7 7 7 7 7 7 7 7 
| V M 113 8 8 8 8*818 8 8 3 8 8 8 84 8£ 8 8 b 6 8 8 8 8 3 8 8 8 8 8 8 
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The core requirements for FUMOKI are 'about AO K bytes and the compilation 
time increases apptoximately with the square of the number of delayed neutron * 
groups-, all other parameters i>eing the same. 

The output of the program is simple. The time, neutron density and 
G matrix can be read out. Also a plot of Che density versus time is 
printed out. 

Problem 1.6.1 \ ' • 

J 

Run the sample problem, i.e. the data shown. 
Problem 1.6.2 

3 

Find the neutron density if you start with 100 neutrons/cm for 

239 

a Pu reactor with all sioc ^delayed, neutron groups. Obtain the solution 
for a 5 second time span and do the same problem with H = 0.01 and (f.l sec. 
Cbmpare the results. Do not write the G matrix. 

Problem U'6.3 y 1 / 

239 

Run FUMOKI for the case for thre^delayed neutron groups for Pu, 

-4 

with a neutron generation time of 10 sec and the reactivity 



p(t) = $0.25 sin 5t. 



Run it for ^a^TotaJ^ime of 1 second, time intervals of 0.01 sec and an initial 
.relative flux of 1.000. Notice that the reactivity at 0.5 sec is $0.17 
and the power has risen to 1.275. 

'If you do 0 the exact same, problem as above for 2 delayed neutron groups 
notice that the* power has become 1.268 at 0.5 sec. / , 



Problem 1.6.4 
• m 

239 



239 & t. 

Given a reactor fueled with Pu and assume that you wish to. use a 

5^delayed-neutron-group model. For a relative power -of .1.000 initially, 

show that after orie>econd the power ia'2.157 if a rod is being removed 

sych that the reactivity inserted is 



' p(t) - P o (l + 5t) 



* All other parameters are as in problem 1.6.3. 

» i • 

Problem K6.5 V * - 

; • N ■ 

Assume a 2 ' 35 U fueled system with a reactivity variation of . \ 

'■' " jT(t) = $0.25 sin 5t. ) 

If the reactor is' to be simulated for three seconds after the reactivity 
is, inserted and..if all parameters, are the jjme as in example 1.6.3 except 
^tor the number of delayed neutron groups then: 

s 

a) use the two delayed neutron model to obtain the- power as a • 
function of time 
and " 

" * b) use the five delayed neutron model to obtain the power* 
variation w,ith time; » 

. How much longer did the computer take to do case v b) than case a)? 

. - ' t 

Answer: About 10 times as long. , On the IBM 370, the times were 11.03 

. sec and .97 sec respectively. Notice also that the power at 

-the end of 3.0 sec is 0.8154 and 0 18161 .for the 2 and 5 group 

f 

cases respectively. ( t * 
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\ % * List of Symbols f or FUMOKI " 

The following syn^ols are listed in^he order of their appearance in 
program FUMOKI . „ 

IFUEL Type of fuel (" 5 U or 23? Pu) (5 or 9) 

> 

NN -Number of delayed neutron groups 



NRO p(t) » p * Type of reactivity to be. inserted. 

p°sin b~t constant reactivity 

p°(l + 5 t) 2»f sine insertion 

0 3 »• ramp insertion 



RO p ~ Reactivity inserted 

.o - 

r 

Bl B2.B3 b-,b 0 ,b 0 Rate or period at which reactivity is 

. 12 3 9 m lnge * te<L } ' 

X(l) X precursor decaj^ constants 

XL A "nAitron generation time 

Time tota J L time tlle reactor transient (s) , 

i is to be simulated 

H At' the 'time step 

jq • t output option to print the G matrix 

- r after 1st iteration(£o print any other 

V " iteration charige statement MKI 2520) 

XN ' n ' initial relative power 

3(X) 'fraction of neutrons for each group- 



5 



BB 
TH 



'total delayed, neutron fraction 
cumulative time used in TPLOT 



11 \ ^ .total number of time steps H fits r 

. . . * - • into TlME * 



r p(t> - | .Reactivity &t any. time 

' , - MM nuffiber of coefficients in the inhour . 

• «•* * ( equation. ; <number of delayed neutron ; ' 

it.- ; • • ; - .': .* rou P s + *> - ' . * ■/ 



H ERJC , .. < ' ' 



., is 



- 42 - 



\W(I). 



M 

WO 

XNCO(I), 
GNXO(I) 

XXNC2 

* „ 

GG.MAD 
MOD (I) 
Z1(I) 

Z2(I) 
AVEBX 

ANHOUR 

GMTRX 
GXN 

POtRT 
TP't^i" 



o ' 



H*(s*ep number) 




coefficients of the inhour polynomial 

degtee* of the inhour polynomial 
(number of delayed neutron groups + 1) 



st eigenvalue of, |A -wlj ■ 0 



column vfector 

logarithm of power 
power 

dimensioned TH used in TPLOT 

precursor density of delayed neutron 
group 1 

reactivity 

subroutine to obtain the average g and 
X depending on how many delayed neutron 
groups are desired 

subroutine to calculate t^e coefficients 
of the inhour equation 

9 

a subroutine which 1 forms the G matrix 



subroutine which multiplies G by 
column vector i.e. G^ 

"A Newton-Raphson iteration 'technique^ 
is used to obtain the roots of the 
inhour equation to obtain 0)^ * 

A subroutine to plot the "vector £ up 
to and including five dependent 
variables. The reactivity, logarithm 
of the power, power and some of the 
delayed neutron precursors are plotted. 



4? 



Flow Chart for FUMOKI 



* Statement Numbers 




9- 2 i?Pu-< 'IFUEL >-5-235U 



1545 



Sub rout 


£ne 




Pu 239 




* 









1650 



Subroutine 
AVEBX 



R01 



I** 



Constant 



Reactivity 





i 



Subroutine 
U 235 



J 






R02 




R03 




sine \ 




Ramp 




variation 




insertion 



r 




t 


1920 


ANHOUR 





Obtains six group g's and 
X's for isotope. * - 

Calculates average g's 'and 
X's depending oh the number 
of delayed groups desired. 



The type of reactivity 
variation desired determines 
which, path to take. 

Routines which determine 

tfte reactivity as functions * 

pf time. f * 



Subroutine which calculates 
.the polynomial form of the 
lnhour equation. 



■ 0. 



7990 



2180 



7880 




) 



Subroutine to determine the v 
] largest root of the inhour 
t equation using New^ton-Raphs.on/ 

method. 



Subroutine which forms Initial 
£ matrix; ^ \^ 



Generate 
Matrix 



"Generates the £ matrix for 
all cases except the first 



Forms the G matrix 



Fofms the product G £ 




Plots the power, etc. 



//WATFIV S?*AKT»GNEGA»PAGC$*du»TIME*i50 



C 
C 
C 
C 
C 
C 

c 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

X 

c 
c 
c 



MAIN PROGRAM 



KINETICS MODULE, 



CODE NAME F M M «U I • * 

OBJECTIVE TO otUAlN IK KINETICS EQUATION WITHOUT FEEDBACK ANO 
SOtVE THE KINETTCS EQUATIONS NUMERICALLY FOR ONE t TO 
•SIX DELAYED N~UTRON^GRdUPS FQR TIME VARYING REACTIVITY 

INSERT 10US ' 
' THE KINETIC ELATION OEKlVED FROM ONE GROUP DIFFUSION 
EQUATIONS • 

FEEDBACK EFFECTS ANO SPATIAL DEPENDENT EFFECTS ARE NOT 
INCLUDED . . «* 



C 
C 
C 
C 

e 
c 

r 
* 

C 

. c 
c 

'c 
c 

t 
c 

.c 
c 
c 
c 
c 
c 

■ c 
c 
c 
c 
c 
c 

, c 

. c 
c 
c 
c 
c 



MK1 
MKi 
MK1 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 



ASSUMPTION- f HE BA^ IC XoUMPTI ON OF THE METHOD IS THAT NEUTRON AHJMK1 
PRECURSOR OFNSITIES BEHAVE EXPONENTIALLY^ WITH A FRE4UENCYMK1 
CHARACTERISTICS OB THE ASYMPTOTIC FREQUENCY CORRESPON- MKI 
OlNli TO THE* REACTIVITY . * * > 

4RITTFN IN SINGLE PRFCISION 



PROGRAM 



GENERAL UESCKl PT ION- UF PARAMETERS 



SYMBOL - 1N/0UT/V DECEPTION 



I FUEL 



NN 



NRC 



Bl 



XL 

R 
A 



WO 
T 

G 

XN* 
C 

TlMF 
H 



IN 
OUT 



IN TYPF I'F FUEL 

5 U-235 

9 = PU-239 
IN* GROUP DELAYED NEUTRON • 

FRCH 1 TO 6 , 
Ity , * %JYPE OF REACTIVITY 
I * \i CONSTANT ^ 

- \ ttU+pm* 
IN V CC\STANT B> .OR .33 

. IN NRQ 
OUT • FRACTION OF DELAYED 
• 'NfJiTROfo 
OUT /?ntu*$GR DECAY 

• < CONSTANT ° * 
IN /NEUTRON bENEPAT IAN 

jf^& , ' 

V V—'rff ACTIVITY 

V' CCtFFlCItNJ OF IwioOft 

FORMULA IN PptY^OMl AL*«^ 
-FOR* . * 

V 1 WFCIP^aCAL' OF ^PERIOD 

V 9 STArtLfc* PLRIOD * u 

V G - MATRIX'* 

V ' VFUTRIJN OENSlfY 

V P*£CUKSP»* DENSITY x 
.N TOTAL TIME 

n TIMC i:>!CfcEM|NT ' * 



- INPUT 

- OUTPUT 



REAL/INT. 



' I 



UNIT 



i/SEC 



1/SEC * 



SEC 



R 


< i/SEC^ 


R 


SE,C . 


R 




R 


N/CM**3 


R 


N/CMM3 


R 


SEC 


R 


SEC 



MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MK 1 
MKI 
MKI 
MKI 
MKI 
MKI 
MK 1 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
MKI 
M^l 
MKI 



20 
30 
40 
50 
60 
70 
80 
90 
100 
1 10 
120 
.130 
140 
150 
160 
170 
180 
190 
200 
210 
2Z0 
230 
240 
250 
260 
270 
280 
290 
300 
310 
320 
330" 
340 
350 
360 
370 
380 
• !390 
400 
' 410' 
420 
430 
440 
450 
460 
470 
480 
490 
500 
510 
520 

53a 

540 
550 
560 
570- 
' 580 
590 
600 
4 ! 



• \ 

VA*IAbLhS 




MKl 
MKl 
MKl 



run PKfs.RA1 IS LIMITED TO DELAYED NEUTRONS FROM THMK1 

TOSS FlSlS^f 1*8 AND PU-239 FOR ANY NUMBER OF MKl 

MKl 

U235.PU239/AVEBX 
K01.R02t A(JO R03 . 



GROUPS FROM GNF TO SIX* 
SUPPORTING ROUTINE 



MKl 
MKl 
MKl 
MKl 

^ MK I 

Sfj^'nt/:... MU-OI. OOJ1CO... -0.10O.. «0U00.. XXNC2U00-K1 

2(6) MKl 

MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 



-OtriNC REACTIVITY AS A FUNCTION OF TIME 



Ri^(KO,T,0L2>^R , i*SIN(ft2*T) 
RJ3<<U\T..d3)*Pa 3 Ml.O*B3*T> 



INPUT PAPAPETfP OAT A 



IIC 



13C 



14C 

15C 

loC 



17C 



1FUCL 



NKG 



XL 

T Mt 

H 

IG 



XN 



1 




TIMEU . 



f- 1 T HCrl 5 OR 9 
^ --FUfcLLD U-235 
HJlLED g-239 
*;PRJP DELAYED NEUTRON , 
TYPE OF REACTIVTTY 
- CONSTANT REACTIVITY 
7 - RtACTIVlTY AS A FUNCTION 
* - Rf*SIN<B2*Tl ' 
3 - REACTIVITY AS A LINEAR FUNCTION OF T 

K = Ra*(l*B?*T) > 
C^STAnT / INITIAL REACTIVITY 
HITHER OP H3 IN NRO 
KEUTrtON GENhRATlCN TIME 
T JTAL TIME USED \ 
UMF INTERVAL 
<t,TPUT SIGNAL F0& G-MATRIX 
1 - PPINT OUT G-MATPIX 
Z DON'T . £ 
INITIAL RELATIVE FLUX^ 



/ 



«v**IMH/T REACTIVITY IN OOLCAK UNIT ******* 



»*F AO 

READ t^> 1 110) IHJFL 
PC AD (^#120) N* 
HEAD (t»f I3C) n OfH'.Bl 
READ (!>.l<fG) XI 
READ (^f)bO) TIMS f H 
ft E A.i> Vj% 16^1 IG 
RE AO ti> 1 170) V 
MR I Tc t6,20l 
• WRITE (6.3C) IFUCttNN 



Thr I NPUT DATA t STEP I P • 26 



617 

620 

630 

640 

650 

660 

670* 

680 

690' 

700 

710 

720/ 
730 

74Q 

750 

760 

770 

780 

790 

800 ' 

8 id 

82C 
830. 
840 
850 
860 



1* <J«J w 

MKip 8.70 i 
MKl '880 



MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 



890 
900 
910 
920 
930 
940 
950 
960 
970 
980 
990 
1000 
1010 
1020 
1030 



MKl 1040 
MKl 105C 
MKl 1060 
MKl 1070 
MKl 1080 
MKl 1090 
MKl 1100 
MKl IUC 
KKl 1120 
MKl 1130 
MKl 1140 
MKl 1150 
MKl 1160 
MKl 1170 
MKl 1180 
MKl 1190 
MKl 1200 



* 4 



51 



1$ 

16 
17 

U 
19 

20- 



r 2I 



22 
•23 

24 

25 

26 
27 
28 
29 
.30 
31 
32 
33 
34 
35 



36 
37 

3a 

& 
40 
<♦! 
42 
43 
44 
45 
46 



*7 

48 
49 

50 
51 

52 
53 
9<* 



55' 
56 



2 J 

30 

40 
50 

oO 

70 

no 

90 
100 

no 

120 
itt 
140 
i5C 
I6u 
L7u 



180 

190 

200 

C 

c 

c> 
c 



210 



*2G 
2 30 



C 
C 

c 

2 40 



wRITc (6,40) XL 

WRITE (-6,50) NKl ,P* V M>U,H1, 

HftlU (6,60) TU'r,H 

WRITE ib,fO) Iii,XN 

WRIT* (6,100) ' 

FuRMAl 1 1HI ,6 Jmm***»* 



f< >o^* ******* ****** t / t IX ?|J6H 

2METEK , / ,2 JH** r ** * >utn*«**nt** f /////J 

f*(JK M A T (2X,liJH TVi>> if mil U-23, I I ,// t2X,?9H NlJMOFP HF 
lUTr in = ,!£,/> 



FQPMAT (^X 

lt>X,»M» ,11, 
FuRiAT t2X 
<NT»»VAL = 



2/H 

22H 



- A7 - 

MKl 
MKl- 
MKl 
MK I 

I . MKl 

*.v4#*«****94***4* KINETICS*' MODULE I **MKl 

INPUT PARAMK1 
MKl 

DEL AYFl Nf MKl 
I* I 



NtUUuij orNLmilUN TIME = , F I 1 . 5 ♦ 2 X , M SE C ) » , / ) .MKl 
IVP»'*;i H ACTIVITY = ,!2t//t?X f 10H Rfl = * ,Fl5.i,MKl 
» = »,Hl«i)t7) , m MKl 

,UH T'.'TAl TlMfc USCO = » F 10 • 3, 2X , f ( SEC > * t 5X , • , TIME IMKl 
',F7.3,» (SEC)*,/) ' MKl 



K'MAT (2Xtl7H OUTPUT )PTION = , I 2 t // t 2X ,25H INITIAL RELATIVE F L UXMK I 
J * ,Fl.3t//> * MKl 

FJK^AI (5X.22H SIX G*!'UP LAMHUA f E ,2X ,6( f 7.4,2X) ,// ) 
FPKMAT C»Xt20H SIX (.KCUP BETA APE ,2X,6(F7. P 5,?X) ,//) 
F'^-IAT (///// f 39H*< s »******* £ Ml) (F INPUT DATA **********,/////) 
Fvi^MAT (L>) ' * 

. FOP MAT (12) 

f Or MAT ( ! IC ,2 U I ;.,:>) ) ' 

FORMAT (FIG.o) , v 

F«>M-AT" ( 2 ( F I J • .> ) ) 
FJMA1 (12) 
F J*- ^AT (F10. 5) 



-CtVPOl- * JURATION AND T YPf. 01 



fuK 



Ml* I 

TH=O.C 
Hl=H ' 

LL^(TlM t *Hl/2.) /HI 
IF (IFUEL-5) l l )JfloO, 
CALL U2J5 (B,X) 
Go TO JOS^ 
CALL POP 59. (o,X ) 
rfkIT? (C,8C) (X( I), 1=1, 6 J 
w*ITE (> ,9o) (> ( i ), 1=1,0) 



•19( J 



-Cr w PtlTf* THt 



AVE KAttf 



iftFT A AND LAMOA F(,P THE 



FUEL 



CALL Avr >X ( { WX,6,,VN'r 
W«ITc (6,210) 

FORMAT (/,5X,27h I HI AV^KACP 

I / .V) 

DO 220 L 3 1 1 NN. 

K*IU (6,233) l. ) ,L,X(L) 
FuP'iAT U X, 2HU( , !2,3r-») = ,Fll, 
00 460 K=1,LL 
G'g TO ( 240,250,.?^C) , krl 



liFTA A NO LAMUA»/»£X<#23H- 



5,5X, 7HLAM{iDA( ,(2,4H) = ,F11.5,/) 



;f t rCtAT UN 1.5.1 



K = h J 
GO T i 



27: 



MKl 
MKl 
MK I 
MKl 
MKl 
MKl 
MKl 

Mki 

MK I 
MKl 
MKl 
NKl 

. * MKl 
MK I 
MK I 
MKl 

\ MK I 

MKl 
'MKl 
HK I 
MKl 
MKl 
MKl 
MKl 
MKl 

PA^TICULAMKl 
MKl 
MK I 
MK I 
, MKl 
-MKl 
MKl 
^Kl 
MK I 
MKl 
MKl 
MK I 
MKl 
MK I 
MKl 
MKl 
MKl 
MKl 



1210 
1220 
1230 
1240 
12S0 
1260 
12TJ 
1200 
1290 
1300 
1310 
I \2Q 
1330 
1340 
1350 
1360 
1370 
1380 
1390 
1400 
l'410 
1420 
1430 
1440 
1450 
1460 
1470 
I4d0. 
1490 
1500 
1510 
1520 
1530 
1540 
15^0 
1560 
1570 
1580 
1590 
1600 
1610 
1620 
1630 
.1 640 
1650 
1660 
1670 
l£80 
1690 
1700 
1710 
1720 
1730 
1740 
1750 
1760 
1770 
1780 
1790 
1800 



ERLC 



57 

58 



59, 



6C 
.61 
62 
63 

64 

65* 

66 

67 

6ft 

69 

ro 

71 



rz 

73 
74 
75 
76 
77 
78 
79 



C 

c fc 

C 

c 

.c 

260 

C 
C 

c 
c 

270 



260 



-see FOUAT ION 



290 
300 

C 
C 
C 
C 



310 
*40 

350 
C 

i - 



R«KO^<wn t TH,Bi) 
GO TO 27J 



•?eC fcOUATIQN l.iu3 



WjRU3(«0,TH,(H) 



MK1 1610 

MKl 1820 

MKi 183,0 

MKl 1840 

MKi 1850 
MKi' 1860 

'MKi 1870 

MKl i860 

*MK1 1890 



CCrfOUlfc THE COEFFICIENT OF INHUU& FORMULA t SEE E QUATMK 1 19pO 



INHOUR FORMULA f / f 4X»35H- 



MKl 1910 
MKl 1920 
MKl 1930 
MKi 1940 
MKl 1950 
MKi i960 
MKl 1970 
MKi 1980 
MKl 1990 
MKi 2000 
MKl 2010 
MKl 2020 
MKi 2030 
MKl 2040 
MKl 2050 
MKl 2060 

LLMPilTl Tie k LjC T Of THl POLYNl «M I At WITH MM COEFHCIEMKI 2070 



ON U3.U WHhkE NN=1 

9 

CALL ANHOUR ( XL iX,tt,*,A,NN) 
MM=NN+2 , 

IF iK.T.Q. 1 ) „WRI T£ (6i230) 
FORMAT (/ f 4X f 3?H ( 01 FF I C 1 1 NT S OF 
! L ,/) 

DO 3D0 !=i,MM 
MW=MM-I*l 
A W ( MW ) = A ( I ) 
IF (K.GT .1 ) Gf 
WRlTfc (6,290) l«AW(>Wl 
F0PJ1AT (/,5Xi<?HA( , 1 2 , 5H ) = 
60NTINU r 
M=NN*l 



,Ml.4) 



AND DfTfcRrtlNC THE LARGEST EIGEN VALUE WO IN STEP 3 



CAjl rul»<T ( A ?(*. Of , M t W , HUUT I t 1 1 R't MM ) 

*Q*'iii I) 

00 310 J=ifM 

V0= AMAX 1 ( W ( J ) iWD 

WO = VO' ' * • 

IF (K.NE.l) GJ TO 370 
WRITE (6,3*00) 

FORMAT (//t^X, 1 THf .LAKGLST t- 1 GEN-VAllUfc 



• »Ell.4,//l 



iCM THH INITIAL VECTOR COLUMN PHI . 

INITIALIZE N fc U 1 R ON DENSITY, T*0 A^O ' PRFCUR SOR DENSITY. 



80 




XNCO( i)=XN 


81 




00 360 I*2tM 


82 




J1=I-1 


83 




8B=»rf+»( Jl)» 


84 




CCJ1 )=B( Jl )/(X( Jl)"*Xl ) 


85 


360 


t XNC0(])=C(J1)*XN> 


86 ^370 


I F ( K .GT , 1 )GG(K)^GNXCU( i) 


87 




GG(K)=XNCO( 1) 


66 




kP=R*B* 


89 




IFIK-1 )330,40vt3fJt 


90 


380 


H=hl 


91 




DO 390 JI=iiM 


92 


390 


XNCOIJl )=GNXCO( Jl) ■ 


93 




GO TO 41) 


94 


400 


H=0.0 




C 






C 






C 


1 .4.15 AND STEP 




C 


TRIX bY. INITIAL' 









f, MAT K I X CORRtSPONOING TO EQUATION 



C 



XKl 2000 
MKl 2090 
MKl 2100 
MKl 2110 
MKl 2120 
MK V 2130 
MKT 2140 
•MKl 2150 
MKl 2160 
MKl 2170 
MKl 2180 
MKl 2190 
MKl 2200 
MKl 2210 
MKl 2220 
MKl 2230 
MKl 2240 
MK 1 j^Wf6 
MKl 2260 
MKl 2270 
MKl 2280 
MKl 2290 
MKl 2300 
MK4*2310 
MKl 2320 
MKl 2330 
MKi 2340' 
MKl 2350 
MKl 2360 
MKl 2370 
MKl .2^80 
MKl 2390 
MKl 2400 



\ 



53 



lERlC,;-; ; ■■• . •• . .• - ' >. 



J 



95 
96 
97 
*98 
'99 

100 
101 
102 



103 
10.4 
105 
106 
107 
103 
109 
110 
111 
112 
113 
114 
115 
116 

117 
118 



C 

410 



^20 
4 30 

440 



%50 



'$60 



433 

*0C ' 



C 
C 

c 
c 
c 



. 49 - 

CML tMTRX (HRiP^X.XL t*UtM,H,G) 
IF (K.NC.2.AN0.IC.NUI) GO TO 430 
H^ITS (o,4fc0) 
DO 420 l*l»« 

WRITE (6,490) tt.M 1 1 1 1 1 MM • *> 
CALL GXN (G»XW.0,^,<;.NX.G0) 
IF (K.LQ. I ) WMTL- C6it*C.) 



v KKl 



MKi 2410 
MK1 ^420 
KKl 2430 
2440 
MKi 2450. 
MKI 2460* 
MKI ?470 
MKi 2480 



O.MAT 7 ,2X :MX ^TI^t^xaHR.TX.ShPO-tR.SX.IZ^X.ZlHuRuUP MKI 2 WO 
LoUAYtb MOTRIA./.IX.hHI » . 1X»0HC SFC I . 5X, 3HC *l . 5X. 7H < - I.7X.10PK1 2500 



MKi 2510 
MKi 2520 
MKI 2530 
MKI 2540 
MKI 2550 
MKI 2560 
MKJ. 2570 
MKi 2580 
MKI 2590 
MKi 2600 
KKl 2610 
MKi 2620 
MKI 2630 
MKI 2640 
iLCUl AT I V. ***MKl 2650 
MKT2660 
MKi 2670 
^ MKi 26«0 

MKI 2690 
MKi 2700 

i( i h ,,n|M IVhO IS CAUULAIIMj tHf AVTWAoc T A ANOMK 1 2710 
<|1h T«" PASK C.F SIX f.WW, . JJKl ?720 

SUPPORT! No MUTIM , Ah I TA 2740 

KKU2760 
MKI 2770 



2H(P*:l.Anve) t//) 
*R1T£ (6,450) kiTliiti IbftiXCCU ) ilMtM) 
FORMAT l3»*X,f 7.3.*X,F 5 .<r . 7< 2X, El 1 .4 M 

IH*TM*H W 
KU)CK|*TH 
UG^GG(K) 
MAI/f-»O=0C> 
11 (K)^GNXC0(2 )]p 
72(K)*P 

XXNC2U) =ALJO( JG) 
WRlTt (6,500) 

CAIL TPL r T (ML0t^Ali»/l , XXNC ^ , t M. ) 
' HiRMAT CVt6Xtl4H Tiu i.-rAT*!X ,/,0X,i5H 
F-UKrtA ! (2X,4(Lll.WX),/l 

FO*>1A*T (//,5X,7C'l< <♦<*•****♦**•******•*«•*** 




119 

120 

121 

12<? 

123 

124 

125 

126 

127 ' 

128 

129 

130 . 

131 

1.32 

133 

134 

135 

136 

137 

T33 

139 



* bUt'fJIJlIlL AVD** |.\,X. J.iU^iUP) 

OIHFNSMN M6). X(<>). UXIM, XX("6), YB(6)» YX<6) 
- 00 13 ^M.N « , i 

XX(K)^X(K) 
U bX(K)=f>(K) 

* GU TO |?5.30t6Jt«*t ili^il^'li f.(». UUP 
20 YC( 4'jFOUP) = ABCTA(hi^M 

Yx( N'jKLMJ" ) = AL AMOA ( »W > t^) 
Go TO 15J 
30 Nl=NGfOUP«-l ^ 
D»J 53 JsUNGROU^ 
YM j) = A8CTA(B,Nl I 
YX( J)=ALAMDA( B, X,M ) 
,hj 4) J 1 = 1 tfU 

\ XI J1)=XX(N2) 
40 »( 
50 LlfrTIMJc 

V GO TO 1^1 * 

60 iO«NGfr°i|i»-l 
J3*2 * 



KK I 
MKI 
MKl 
MKl 
MK I 

- MK I 
MKi 
MKl 
MK I 
MKl 
MK I 

. MKl 
M 1 
MK 1 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKi 



2 780 

2 700 

2«00 

2H\0 

2820 

2830 

2840 

2850 

?<)60" 

2870 

2880 

?a90 

2900 

2910 

2920 

2930 

2940 

2950 

2960 

2970 

29B0 



1 
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ft' 







DO 90 JL»1*NG'*0UP 


i"4o . 




141 


, * 


YB(Jt)rA6eTA(BtN3) 


"142 




YXJ J 1 ) = AL AMDA ( B ♦ X*N3 ) 


♦ * 
'143 




If .Ul.EG.NGROUP* GO 


144 




00 70 J2=ltN3 


145 ' 




J3^J3 + 1 


146 




B( J2)=6X( J3)., 


147 


70 ■ 


Xi J?)=XX< J3) 


148 


80 


CUNT I N\) 6. 


149 




GO TO 15o 


150 


90 


N4 = NGP0U<>-2 


151 




J3=2 


152 




00 UC Jt = i»N4 


163 




,YB(Jl) = A3ETA( Bt N4 ) 


154 




YX(Jl)=ALAMDA(BtX.N41 


155 




00 100 J2=1.N4 


156 




J3 = J3 + l 


157 




B( J2)=BX(J3) 


158 


100 


XU2)=XX(J3) " 


159 


1 10' 


CONT 1NUE 


160 




00 120 N')=3»NGROJP , 


161 




NI=N0+2 


162 




YB(NO)=BX(NI ) ^ 
YX(NO)=HX(Nl) ' X 


163 


120 


164 




GU TO 150 


165 


' 1-30* 


N5=NGPDUlM 


166. 




Yii(l)=AB^TA(B,.?) 


167 




YX(l )=AL\MDA(B,X,<>) 


168 




00 140 Jl*?tNG*CUP 


lo9 


t 


J2*JL+l 


1 70 




Yb( JlJ-rX( J2) 


171 


140 


YX(J1)=XX( J2> 


172* 




GO T0U5 ) < 


173 




JO '160 J2=l»\GkCUP 


174 




» tf< J2)=YBU2> 


175 


16C 


X(J2)=YX( J2> 


176 




KBTORN 


177 







50 



S 



c 
c 
c 



Pim.TION ALAKDA IS AVERAGING THE VALUE OF LAMriDA 

SEe EOUATlfN 1.3.5 



1 TH 
179 
180 
181 
182 
18* 
184 
185 
186 
187 



10 



FUhUTUN ALAMO 
0IHENS401 a v tr*l 
B, 1 = 0.0 
8L=0.0 • 
00 10 ! = UN 
Bl=3l+B( 1 ) 
BL=BL*8( I )/X( I ) 
AL*AMDA = B 1/BL 
PETURN 
END 




f-ONCTIOM Af>£ T A IS AVERAGING BETA 

SfE EO.UAT I GN 1.2.3 



flKl 


2990 




3OO0 


MKl 


3010 


MKf 3020 


MKl 


3030 


MKl 


3040 


' MKl 


3050 


MKl 


30<jO 


MKl 


3070 


MKl 


3080 


MKl 


3090 


MKl 


3100 , 


-MKl 


3 HO 


MKl 


f 3120 . 


MKl 


3130 


MKl 3140 


MKl 


3150 


MKl 


3160 


MKl 


3170 


MKl 


3180 


MKl 


3190 


MKl 


3200 


MKl 


3210 


MKl 


3220 


MKl 


3230 ■ 


MKl 


3240 


MKl 


3250 


MKl 


3260 


MK I 


3270 


MKl 


3200 


MKl 


3290 


MK I 


3300 


MKl 


3 310 


MKl 


3320 


MKl 


3330 


MKl 


3 340 


MKl 


3350 


MKl 


3360, 


MKl 


3370 


MKl 


3380 


MKl 


3390 


MKl 


3400 


MKl 


3410 


MKl 


3420 


MKl 


« « 
3430 


MKl 


3440 


MKl 


34'50 


MKl 


3460 


MKl 


3470 


MKl 


3480 


MKl 


3490 


MKl 


'3500 


MKl 


3510 


MKl 


3520 


MKl 


353tf 


MKl 


3540 


MKl 


3550 


MKl 


3560 



} 



&ERIC 



55 
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MK1 3570 
MK1 3580 



PUNCH ,H 
DIMFNSIO 
A8fcTA-0 • 
DO 10 1= 
Aht T A=A6 
RETU^ 
fcMD 



♦ ' MK I 3590 
MK1 3600 
MK1 3610 

* MKI 3620 
MKI 3630 
MK1 3640 
HK1 365,0 

*« MKI -3660 
MK1 -3670 

FUR SIX GRGU* OcLAYtli NCUTRON FROM THfcRMAL FISSION . . MM 3690 
> MKT 6710 




X(6)=3.J>1 

RETURN 

END 



MKL 3720 
MK1 3730 
/ MKW3740 
MK1 3750 
MKI 3 760 
MK I 3770 
MK1 3780 
MK 1 3 790 

\ MKi 3iT0(r 

MK 1 3&10 

' " • MK 1 3620 

MKI 3830 
MKI 3840., 
MKi 385(1 ' 
MKI 3860 
MKI 3870 
MKi 3080 
. MKI 3S90 

^hkilUTI*'. PU-?39 IS PRUVjniM, DATA F-k BF.TA AMU LAMHK1 3900 

OA Sit r,«OU* OJbLAYFP NFUTRON FPC* THIRMAL FISSION . MKI 3910 

MKi 3930 



SUft^.JUl INfc U23 5 (B.X> 
UjyjfcMSlOf/ H(o)t X(6) 

£>(2>=. 00141 » 
• fr<3)=. 00127 
rt(4) = .0O<>55 

i3 ( o) =.00027 
X<1)=.0124 
XK2)=.0305 
X(3)=.lll 
X(4)=.>01 
X(5)=i. L\ 



SUtfKOUTINF PU2^ 

CUMfNlSl^N tt<6>. 

3Cl)=,OCC094t>0 

»( 2)=. 0 5030460 

tf(3)=. 00056970 

b(4)=.000Q8J20 

3(5)^=. 00023220 

p{6 J=.000118dO 

XU J-.0128 

X(2l=*0301 

X(3>=.124 

X( 4) = .*325 * 

X(5M1.12 

X(6)=2.69 

RETUPN 

lNO 



X(o> 



MKI 
MKI 
MKI 
MKI 
MKI 
.MKI 
MKI 
MKi 
MKI 
MKi 
Ml$l 
MKI 
. MKI 
MKI 
<MK 1 
• MKI 
MKi 



-mo 

3950 
396C 
}970 
3980 
399^1 
4000 
4010 
4023 
4030 
4040 
4050' 
4060 
4070 
4080 
4090 
4100 



56 



- 52 - > 

* \ MKl 4110 

THb MAIN PRCXtRAM IS TRYING TO CALCULATE THF C06FFICIMK1 4120 

ENTS Of INHCUR FU8MULA FOR ANY JGROUP B6TWEEN 1 TO 6 . MKl 4130 

*i MK I, 4 140 

* * MKl 4150 



SUBROUflNL ANtHU" ( XL t X » 3 , K , A , NN) 
D I MENS I W X(NN)t MNN)i A(8) 
B»=0.0 

DO 10 l^itNN 

B6 = 8rt + IM I ) ^ * 

J R0=R*BB 

JJ=NN+2 

DO 20 J=1,JJ 

A(J)=0.0 c 

Sf E L1V.1 3 Fi.f. NN=1 
GO TO ( 30»40,^0f6C,70,a0) , NN ' 
CALL Gi (XLtX,rt»R(.»A); 
GO TO 90 • 

CALL 02 (XLtX*U»Kf t A) * ' 
GO TO 93 . 

CALL G.J (XLtXttitKl t A ) 
Gti TO 90 

CALL G4 lxTtXtH««('t A) •* * 

GU TO 90 

CALL G5 (XLtXtbtRtSA) 
S>U T J 90* 

CALL Go (XLfXttlt*('f A) 
GO TO 
RETURN 
cNO 



s 



MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
, MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
— ^ MKl 
MKl 
" MKl 
MKl 
^Ki 
-MKl 
MKl 
MKl 
MKl 
" MKl 
MKl 



MKl 

SIM^UJINfc 61 »G?»G?tG4»G5»G 9 t IS 01 RECTL* C ALCULAT I MK 1 

COEFFICIENTS OF INHUUR FORMULA IN POLYNOMIAL FORM IN^ 

*t£e ORDER CF THE smallest DEGPEF TO THP LARGEST • 
|Ht " FORM OUTPUT IS A1*W**.N'+ A2*W**(N-1) + A3*W**(N-2) 

+ ♦ ,an 

SCreOy*»TICN U3.1b FOP THE GENERAL INHOUR FORMULA 

SUPPOrTiNi, PCUTINt A62tA52tA4?>A32,A22 AND 
. " XX,XY*XZtXUt xv,a£ »•> 



MKl 
MKT 
MKl 
MKl 

Vki 

MKl 
MKl 
MKl 
tfKl 
MKl 



4160 

4170 

4180.^ 

4190 

4240 

4210 

4220 

4230 

4240 

4250 * 

4260 

42T0 

4280 

4290 

4300 

4310 , 

4320 

4330 

4340 

4350* 

4360 

4370 

4380 

4390 

A400 

4410 

4420 

4430 

4440 

4450 

4460 

4470 

4480 

4490 

4500 

45J>0 

4520' 

4530 



.SUbROUTINE iii .(^LtX^OtRPt^) 
"L I MENS! ON XH)t H< 1 ) i *<3) 

tt(2)=Vi 1)*XL*X< i)-fC 

H( l)=-KiVX( 1) , . 

P6TURN 

END * 1 



MKl 4540 
MK l " 4 550- 
MK1 4560 
MKl 4570 
MKl 4580 
MKl 4590 
MKl 4600 



SURRUUUNfc G2 (XLtX«BtRCtW) 
DIMENSION X(2)t tSi/U W(4) 
W(4>=*L 

h(3) = XX(8t2)fX'L*XX( X,?)-RO 



MKl 4610 
^MKl 4620 
MKl 4630 
MKl 4640 



57 




2*>2 
263 
264 
265 



266 
267 
268 
269 
■27Q 
271 
272 
273 
274 



£ 275 
276 
277 
278 
279 
280 
281 
282 
283 
284 



■ ^ - 53 - 

wm~Xl*AAtX» 3) +A62<B,X,,2)-R0*XX<X,2> 
W(U*-KO*AA<X*2 ) 
RETURN 

f NO < 



SUBROUTINE G3. ( XL , X, b , *0t W i 
DIMENSION X(3>» IM3>, WJ5) 
.W(5)-XL 

W(4)=XX(B.3 ) + XL 'XX(X f 31-KU 
W( 3)=XL*XY( X, 3)*A6?<HtXf J)-POtfCX( X,3) 
> W(2) = XL*AA(X, 3)+Ai>?{htX,3)-RC*XY(X,3) 
W( i)=-RU*AA{X,3) . 
RETURN 
END 



SUBROUTINE 04 (XLfX,ti,PCfW) 
Dl MENS iCN X<4). 6(4), W(6) 
W(6)-Xl 

Vi( 5)=XX(B t 4 )4-XL«"XX( X, 4) -KG 

W(4)=XL*XY(X,4) ♦A62<fi,X,4)-R0*XX(X,4) 

W( 3)=Xt>XZ( X, 4) *;\52<n, X, 4 ) -RO*XY( X,4 ) 

WI2)=XL*AA(X, 4) ♦ A42,(ti*X.4)-RU*XZ(X.4) 

W( 1 )=-R')*AA<Xt4) 

RETURN 

END 



MKl 4650 
MK1 4660 
MKl 4670 
MKl 4680 



MKl 4^90 
MKl 470? 
MKl 4710 
MKl 4720 
MKl 4730 
MKl 4740 
MKl 4750 
MKl 4760 



MKL 4770 



Kl X?90 
K/t 4800 
I 4810 
MKl 4820 
MKl 4830 
MKl 4840 
MKl 435 
MKl W60 
MjpK4870 





285 
286 
237 
288 
289 
290 
291 
292 
293 
294 
295 



SJtMjUTlML 05 (XLfX»!W<*Cffc) 
DI^NSI )N X(5)i r>J5), «(7) 
W(7)=XL 

rf(6i=XX( i,5) + XL*XX( X,3f-KU 
W(5)=XI*XYIX, 5) ♦Ao2(B.X.5l-RC*XX(X.5> 
W(4)=XL*X/( X, 5) ♦A52(fr,X,5)-P0*XY(X,5> 
M 3)-XL*XU(X,5) + A<*?(tl,X,5)-Ru*XZ(X,5) 
W(2)=Xf*AA(X,5) +A32(B r X,5h-RC*XU(X,5) 
W( 1 )=-R»)*AA(Xf 5 v ) 
RETURN / 
END 



MKl 
MKl 
MK 1 
MKl 
MKl 
MKl 
MKl 
MKl 
' MKl 
MKl 
MKl 



4880 
4990 
4900 
4910 
4920^ 
4930 
4940 
4950 
4960 
4970 
4980 



296 
297 
298 
299 
300 
3.01 
302 
303 
304 
305 
306 
307 



iU6RJ0TlN;E 06 (XLtXi8»Kt ,1^1 
DIMENSION £<o*\ l J(())i-«(BI 
w(8)=XL . 

W ( 7 ) - X X ( « t 6 ) t Xt>*< x I X .t o I - A C 
W(6)=XL*X.Y(X,6V + W(^,X',6i-Pt*XX(X t 6i 
W(5)=XL*XZ(X,6) +Ab^(r,X,6)-RL*XY( X,6i 
W(4) -XL*XU( X,6) ♦A4/(fi.X. 6)-Ru*XZ( X, 6) 
W(3 )=XL*XV<X»6) ♦ A3?tB»X, 6)-fa*XUt X,6) 
W( 2)=XL 3 >AA(X,6H-A^2itt,X,6)-RU*XV(Xt6) 
W( 1 ) = -*0*AA(X,6) 
E TURN 

END * 



-fUrjfTH.'N A62 fUNCTION AS : 

A62 s» f l*(X*+X3+X4+»««) ♦ 02* ( Bl +H3 + R4 + . . . . ) ♦ 
B2*(Bl + lU+ft*+*J4+«"> * IN THAT C0M8I NAT I 0H% 



MKl 
MKl 

b'K 1 

MKl 
MK i 
MKl 
MKl 
MKl 
MK 1 
MKl 
MKl 
MKl 
MKl 
MK 1 
MKl 
MK 1 
MK 1 
MKl 



4990 
5000 
5010 
5020 
5330 
5040 
5050 
5060 
5070 
508.0 
5090 
5100 
5110 
5120 
5130 
.5140 
5150 
5160 



308 

309 

310 

311 

312 

313 

314 

315 

316. 

3A7 

318 



10 



So 



c 
c 

« 

c 
c 
' c 
c 
c 



- 54 - 



FUNCTION A62 (BtX.M) 

0IMENS10N X(M), AXl(6) t 

A62=0*0 

00 20, 1 = 1, M 

00' 10 J = l.f 

AXK J)=X( J) 

Wf>><B(I)*XX<AXl,M) 
A6?=A&?*W( I ) 
RETURN- 
END ' 



^ A5? IS ALMOST THfc SAME 

.'hXCcPT THAT XI * YIAY2 WHIC 
EOUAl :N»: IM HWS1 I HAP AC T 
S MAI LE k T'IA.\ Nt XI 



THING AS A62 IN 



MK1 


5170, 


MK1 


5180 


MKI 


5190 


MK1 


5200 


MK1 


5210 


MKI 


5220 


j MK1 


5230 


MKI 


5240 - 


MK1 


5250 


MK1 


5260 


^MK1_5270* 


MK1 


52 80 


' MK1 


5?90 


. ^ MK1 


5300 . 


IT'S FUNCTIUNMK1 


5310 


NEVER *i E MK1 


5320 


At WAY'S MK1 


' 533ft 


MK1 


5340 


MKI 


5350 


MKI 


5360 



319 FUNCriUN A52 TB f X t * ) r * * MK 1. 5370 

320 OlMSNSlrtN X(M>. AX2(<>>, rtff % )t WU) 4 MKI 5380 /* 

321 • A52=0.0 , * MK1 5390 7 

322 00 20 l = ltM, n MK1 5400 

323 00 10 J-1,M *M 5410 

324 10 AX2(JJ=X(J) , MKI 5420 

325 AX2(I)-0.0 . MKI 5430 
326^ W( I) = *< I)*XY(AX2«M) M*l 5440 

327 20 A52*A52*W(I) * ^ MKI 5450 

328 RETURN MKI 5460 

329 END , MKI 5470 

C * MKI 5480 " 

r MK 1*549*0 

c , u \ u AS Ao; EXUTT ITS AUDITION OF THRFE ChAKACTEK ,MK1 5500 

0 MKI 5510 

- C * 'MKI 552Q 

330 F,UNCT Ii)*M A^2 (RiXi"'J MKI 5530 

331 # DtMFNSlOM X ( M ) » AX3 ( 0 ) • ( H ) • w ( 6 J • MKI 5540 

332 A42=0,0 ^ . MKI 5550. 

333 DO 2.0 I-L, H * MKI 5560 

334 < w !>0 10 J*lfM * MKI "5570 

335 lH» AX3(J> = X(J) , . MKI 5580 

336 AX3U) = 0.0 " " ^ ,.*K1 5590 

337 W ( I ) =B ( I ) **XZ ( AX 3 t M ) , 4 ^\ MKI ^600 

338 20 A4£»A{2+WII J o . * \ MKI 5610 

339 RETURN . ; » \ KK 1 5620 

340 END \ MKI 5630 * 
C * * ' ' \ MKI 5640 

0 " * V MKI 5650 • 

c ^ ^r.Q AS, A62 EXCEPT 'CHARAUfcP A00ITI0N IS FOUP . * 'MKI 5660 

C , '•',•» ' * MKI 5670 

, t C 1 , "Ki 5680 



\ 




FJUN0J.K1N A32 <h,X,M) • ' , 
OJ-MENSION X(M)» AX4(6»» U(M),<Wl6) 
A32»0.0 

00 20 l=lt& . ' 

00 10, J* If M 
AX4U)=*XU) 
AX4(I)-0.0 • 

wtn=B«n*xu(AX4,M) 

A3^A32*WU) 

RETURN - ' • 



— wF.*) AS A62' FXCt'PT AOOITICN CHARACTFR.IS FIVE 



MKl 

Ml 

MKl 

MK1 

MK L ' 

MKl 

MKl 

MKL 

MKl 

MKl 

MK I 

MKl, 

MKl 

MKl 

MKl 



5700 ' 
5710 
J>720 
5730 
574Q 
5750 
5 760 
5770 
5780 
5790 
5800 
'5810 
5320 
5330 
5840 



352 

35<r 

355 
356 
357 
358 
359 
360 
361 
362 



10 



20 



C 

5 
c 

c, 

c 



FUNCTION A22 <i«tX,M) 
DIMENSION X(M)t AX5(6) t l>(M) f W(6) 
A22=0.0 ' , t 

DO 20 I-ltM . 
00 10 JM,K 
•AX5U) = X<J) 
A*5(I)=0.0 
*( IKsBf I)*XVf AX&.M 
A22=A22+W( I ) 
R£ TURN 



X-X IS CGMPUT I NG .THE fflTAL Uf- X 

XX - m+x2+x3+x4 + .7..* +XN) 1 



MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MK l f 
MKl 
MK V 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 



5850 
5860 
5870 
5880 
5890 
5900 
5910 
5920 
5930 
5943 
5950 
5960 
5970 
5980 
5990 
6000 
6010 



363 , 

364 

*365 

366 

36 7 

368 

369 



It 



FUNCT I UN XX (XtM 

CIMPNS10*J X<M) 

XXO.C 

00 14 l*l#H 

XX = XX + XU) * 

ReTtJftN 

END 



„__„ XY IS CP1PU1I.NI. THf AOOITIlA CF TWC PKWHJf^S 

XY = <Xl*X2 + Xi*X3+Xl<*44- +X2*X3*X2*X4+. . . ) 

THE FIPST POSTCRIPT 



;?feV«AYS SMALLel^ THAN THE: 



MKl 6020 
MKl 6030 
MKl 6040 
MKl 6050 
MKl 6060 
MKl 6070 
Mtfl 6080 
MKl 6-090 
MKl 6100 
MKl £110 
MKl 6120 
MKl 6130 
MKl 6140 
MKl 6150 



370 

371 

372 

373 

374' 

375 

376 



PUNCT ! "N XY ixt«> 
DIMENSION X<M> 

xr^o.o 

00 20 I-ltH 
OC 20 J^ltM" 
A = XU) i 
^ = X(J> 



I' 



MKl 6160 

MKl 6lft) 

MKl 6180 

MKl 6190 

MKl 6200 

MKl 6210 

MKl 6220 



J - 



it 



60 















1 


• 




: 377 




IF i I-JJ 10i20,?0 


378 


10 


Y -A*B 


379- 




XY=XY+Y » - 


t 380 




CONTINUE - 


, 381 


• 


RETURN 


382 




END 




C 






c 






c 






c 






c 
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Ori) A S XY FXCtPT IT IS A PRODUCT OF THREE 











383 




FUNCTION *XZ ( X,M) 


384 




DIMENSION 




385 




XZ=0.0 




386 




00 30 1=1 


f M 


387 




00 30 J=I 




388 




DO 30 K*l 


tM . 


389 




A=X,( I ) 




390 




B=X('j> 




391 




C=X(K) 




392 




IF <I-J> 


10.30*. 30 


393 


10 


IF (J-K) 


20t30.30 


394 * 


20 


Y=A*B*C 




395 




XZ=XZ+Y 




396* 


30 


CONTINUE 




397 




RETURN 




398 


C 

C 


END 






*.c 
c 

c 


— ~— — - 





K) ^ XY ^XCtPT FOR PRODUCT OF FOUR ..^ 



399 

400 

401 

402 

403 

404 

40t> 

406 1 

407 

408 

409 

410 

411 

412 

413 

414 

415 

416 

417 



r 



10 
20 

30 ' 

40 



C 

c 
c 
c 

c 



FUNCTION XU lXtH> 
DIMENSION X(M) 
XU=0,0- 
DO 40 1 = 1. M 
DO 40 J«l?fc *\ 
DO 40 K = lt* 
DO 40 I =l.M 
A = X( f > 
B*XU) 
C=X(K) 
D=X<U 

If ( I -J) I0.tt«0«40 
IF (J-K) 20t40i43 
IF (K.-U 30V40.40 
Y=A*a*f*D 
XU=XU+Y 

CONTINUE * v , 
RETURN 
5N0 ' - 






MKI* A? "TO * * 




MK 1 A740 








Ml/ 1 A? f»0 • 




MK 1 70 * * 




« MK1 6280 




MKI 6290 




MK1 6300 


OF THRcc • 


MK 1 fc^t 0 




MKI 632P * ? * 




MK 1 6330 




nr\ i* kj j ~ *j 




MK1 6350 




MKI 6360 


* 


MK 1 f*370 




uvi AO Q (\ 
MM O jOu 




MK 1 A^qn 




Pif\i Otww 




MK 1 A41 fl 




MK 1 A4?0 


> 


MK 1 r*4^0 




MK1 A 440 


g 

r 


111/ 1 tACA « 




MK 1 fi4f»0 




MK 1 f»470 




MK I 6480 




MK1 649*0 




MK 1 A500 




MK 1 6510 


FOUR • 


MK 1 A S ? 0 ^ 




MK I 6530 - 




MK 1 f\ S40 




MK f "ft'j'jO 




MIC 1 fy^fxQ * 




MK 1 T»S70 




MK 1 A^AO 




' MK 1 A'iQO 




nf\i oouu 




* MK 1 AA1 0 




mk i f\(t?a 




MK 1 AAlf) 




MK 1 AA40 

PIN I OOtv > 




Mki A^SO 




nf\ 1 OOuu 




MK 1 AA7n 




MKI AAR0 




MKI AA90 




MKI fc700 




♦ MK 1 A710 




MKI 672 0 




MKI 6730 r 




MKI 6740 




MKI 6750 



'DOIT ION OF MULTIPLICATION OF FIVEMK1 6760 

tfKl 6770 
MKI 6780 



1 ERJC 



- 57 - 



U16 
j 419 
j 420 
i 421 
422 
423 
424 
425 
426 
427 
> 428 
' 429 
43.0 
431 
432 
433 
434 
435 
436 
%37 
4?* 
439 
440 
441 



L0 
20 
30 
40 

50 

60 



X 
C 



FUNCTION XV UUM) 

DIMENSION. X(M) 

XV*0,0 

00 50 1=1. M 

oo 50 » ; 

DO 50 K=l,fW 
DO 50 IMfMfr 1 
DO 5& N=l 
A*XU) 
B=X( J) 
C=X(K) 
D=X<L) 
E=X<N) 
I F (I-J) 
IF (J-K) 
IF 

ir 

Y=A*3*C*D*l 
XV=XV*Y 
^CONTINUE 
-WRITE (fa* 60)* XV 
. FORMAT (2X, 7H* XV 
RETURN 
END 



20*50,50 
JK-L) 30,50, i>0 
(l-NI 40,50, 5-? 



' J 



,F19.10) 



AA 



-AA IS COMPUTING MULTIPI ICATION OF N CHARACTER 
X1*X£*X3*X4* *XN 



MK1 

MKl 

MK 1 

MKl* 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

'MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 



6790 

6800 - 

6810 

6820 

£830 

6840 

685>0 

6860 

6870 

6880 

6090 

6900 

6910 

6920 

6930 

6940 

6950 

6960 

6970 

6980 

6990 x 

7000* 

7010 

7020 

703.0 

7040 

7050 

7060 

7070 

7080 



442 
443 
. 444 
:44 5 
446 
447 
448 



10 



C 

c 
c 
c 
c 
c 
c 
c 
c 
c 



FUNCT UN A A t>X 
D I MENS I l V\ X-(M) 
AA = X( 1) 
DO 10 
AA=AA* 
RETURN 
END 



I=2.M» 

* LP 



0 



j: ft* ft k X N MATftlx'bY HANSEN METHOO ., & - "ATT* IX . 

The FI^ST COLUMN , IMF. FIRST ROW AND THE DIAGONAL * 0, 
THE RtSr Of- THFM = 0 . . . 

SfcE L0UAT1CN L.4.15 



SUPPPFTIkO ^UUTINt- NONf* • 



1> 



MKl 7090 
MKl >7 100 
MKl 7110 
MKl 7120 
MKl 7130 
MKl 7140 ' 
MKl 7150 
MKl 7160 
KKl 7170 
MKl 7L80 
MKl 7190 
MKl 7200 
MKl 7210 
MKl 7220 
MKl 7230 
MKl '7240 
MKl- 7250 



449 
450 
45f 
A52 
453 
454 
' 455 
456 
457 



SUtMOUJIKfc GMUX C r , B ,X , XI , WU , N, »> U > ^ ^ • 

DIMENSION b(6)t X(6), btN.N* 
O0(AX,»U = LXP(-AX^H) ' 

Rh(H,BX)=(R-8X)/XL ' , 

NN=N-1 * . 
DO 10 1=1. NN 



MKl 7260 
MKl 7270 " 
MKl 7280 
MKl 7290 
MKl 7 300 
MKl 7310, 
MKl 7320 
MKl 7^30 
YK1- 7340 



IERIC 



62 
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458 
. 459 
1 460 

: 46i 

1 462 
. 463 
\ 464 
.f465 
466 
, 467 
468 
469 
. 470 
471 
472 
, 473 
474 
* 475 
476 
477 
478 
. 479^ 
480* 
481 
482 
483 
484 
465 
486 



48 7 



10 



20 
30 

40 
50 

60 
70 

80 



90 
100 

i 

110 
120 

130 



140 

f 
C 

150 



488 




489 




* 


C 




c 




f 




c 


\ 






c 




c 


V 




490 


* 


491 




492 




493 




494 




495 




49fe 


20 


497 




498* 






c 




c 




c 




c 




c 




c 


* 


c 




c 



00^150 J*l,N 
\D0 LfO JJ=1»N 

RBL*RB<.R,8B) 
IF <J r JJ> 20,60.20 
IF. U-l) 30.50,30 
Jl=JJ+l 

IF (JJ-N) 50,40,50 
IF (J-N) 120,80.120* * 
AX=X(J1-U ' 
GU TO 90 

IF IJ-IJ 80,70,80 
GU,l) = EXP<RBtP,BB)*H) 
GO TO 140 
AX=X ( J- t ) 
G< J,JJ)=GD(AX,H) 
GO TO 140 " 
If { J-l) 110* iOC.UC 
G< l.fJ^)=GH( RBLtWO.H, AX) 
GO T0140 

IF (JJ-1) 123,133,120 
GfJ,JJ)=0.0 

GO TO 1 40 ' , 

AX=X( J-l ) 
BX=B(J-l> ' 

G( Jt 1)=GVIXL» Wp,H,AX,flX) 
* CONTINUE 

.WRITE (6,144) (G( I , I I ) , ! 1/1 • N ) 
144 FORMAT ( > Z ( £ 10.3, 2X0, ) 

CONTINUE 
RETURN ' 
END 



MULTIPLY THl G 

TO GET TMf: NEXT UNI , 
SEE E0UA7 I r,N l.^,J. ; 



• MATRIX WI-T^ THE 



MK 
MK 
M»Si 
MKll 
MKjl 
MK1 
MKl 
t M»f 1 

mi 
i 

^ 

V «MK 1 
J$K1 
PlKl 
MKl 

^ /MK 1 

/ MK 1 
/ MK 1 
/ MK 1. 
/ MK 1 
/ MK I 
/ <MK1 
MK1 
MK 1 
MKl 
MK I 
MKl 
MKl 
MK 1 

/ MK 1 

MK I 

/ MK I 

MK 1 
MKl 

/ MK 1 

L k MK 1 

Initial vector colmki 

MKl 
MKl o 
MKl° 
MK 1 



SUBROUTINE GXN tG , AN , N , GX-J 

DIMENSION G(N ,N ) , XN(NJ, GX(N) 

DO 25 1=1, N 

GG=0.O' 

DO 10 J=l,M ' K 

CG-GG+G ( I , J)*X"4( J ) 

GX(t)*GG < 

RETURN „ v * , - 

END ] - 



SUB&jUTINE POLHT 



7350 
7360 
7370 
7380 
7390 
7400 
7410 
'7420 
7430 
7^40 
7450 
7460 
7470 
7430 
7490 
7500 
751-0 
7520 
7530 
7540 
7550. 
7560 
7570 
7580 
7590 
7600 
7610 
7620 
7630 
7640 
7650 
7660 
7670 
7680 
7690 
7700 
7710 
7720 
7730 
7740 
7750 



1 



PURPOSE , ' 

COMPUTES THF- RcAl AND CIJMALEX ROOTS OF A REAL POLYNOMIAL 



MKl 

MKl 

MKl 

MKl 

M,Kl 

MKl 

MKl: 

MKl 

MKl 

MKl 

MKl 

MKl- 

MKl 

MKl 

MKl 

MKl 

MKl 



7760 
7770 
7780 
7790 
7800 
7810 
7820 
7830 
7^40 
7850 
7360 
7870 
7880 
7890 
7900 
7910. 
7920* 



63 
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USA C4LL PQLRTIXCi^fCoSpiMtHOOTP.ROOTItlFK.Mll . 

ntsr^rPTijN of paramfters r -' 

XCOF -VECTlA 01 A*l COEFFICIENTS OF THF POLYNOMIAL 

3R0£«n FROM SMALLEST TO L ARGF $T POWER ' _ 
COr -WORKING VICTOR OF LrNGTH M + l 

M -ORDER-Uf POLYNOMIAL. 1J , inTC 
- ROOTk-RfcSULTANT VICTOR JF LENGTH M CONTAINING PEAL ROOTS 
OF TH£ POLYNOMIAL 
RPOn-PCSULlAM VKTOR I'F LENGTH M CONTAINING THE 

CORKSPr,NnUlG IMAGINARY ROOTS OF THE POLYNOMIAL 
IFR -E**CP COOC ViHCPL . , , 

I Ek=^ M CRRC** 
I PR = l K LESS THAN ONE 
* IER = 2 M GRf AT E P THAN 36 ✓ , 

re^ = 3 UNAHLC TC DETERMINE ROOT WITH 500 I TER AT I JNS 

ON 5 STARTING VALUES 
IfcA*4 "HIGH ORDER COEFFICIENTS ^ERO 
Ml * -NUMBER UP COEFFICIENT » M+l 

R£ MARKS - rrc 

LIMITED TO 36TH uROtR POLYNOMIAL OR LESS. 
PLIATINC POINT uVERFLOW MAY OCCUR FOR HIGH ORDER 
PUYNUHIAIS .UJT Hi LL NOT AFFECT THE ACCURACY OF T Hi 
RESULTS . • 

SUBROUTINLS A'iL FUNCTION SUBPROGRAMS PIOUIRED 
, N )NE \ * , 

,,:1 nKt0N-kAPHS-1Ii ITERATIVE TECHNIQUE • TMl FINAL ITFfrATlWS 
' ON EACH Ht.OT.ARF PtkFOPHEO USING THE IRlGItiAL POLYNOMIAL 
RATHER THAN THE RlDUCEb POLYNOMIAL TO AVOID &CCUMMLAT tO 
eRRORS IN THE REDUCED Pf L YNOMJ AL . 



MK I 
MKl 
MK i 
MKi 
MKl 
MK \ 
MKl* 
MK 1 
MKi 
MK 1 
MK 1 
MK 1 
MKl 
MK 1 
MKl 
MK 1 
MK 1 
MKl 
MKl 
^Kl 
^MKl 

\k I 

MKl 
MKl 
MKl 
MKl 
MKl 
* MKl 
?^K1 
-MK 1 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 
MKl 



7930 
7940 
7050 
7960 
7970 
7980 
7990 
8000 

ana 

8320 
'3030 
8040 

8050 

8060 

8070 t 

8080 

8090 

8100 

8110 

8120 

8130 

8140 

8150 

3160 

8170 

8180 

8190 

8200 

8210 

8220 

8230 

3240 

8250 

8260 

8270 

8280 

8290 



SUBROUTINE POL>J < X CUP , LUf . V , MIOT R , R'JU I . I P , ,M ) 
UI^ENSION XCOF(Ml), CuMMl), RUOTR(M), ^Jl(M) 

00U8L E PRECIS* XII. Yf. X. Y . XPk . YPP •UX.UY. V. YT f XT.U .XT2 ,YT2 . SjJMS'J. 
1 DX , DY » TEMP i ALPHA f CAMS 

IF A^OOUBLE PRECISION VbPSION OF THIS ROUTINE IS DfcSIJbO, THE 
J ill COLUMN I SHOULD BE Rf MOVE D FROM THE DOUBLE PRECISION 
STATEMENT WHICH FOLLOWS. ' ' 

DOUBLE ^RtClSION XCCF iU'F ,*COTP .ROCTI 

THE T MUST' ALSO Bl RCMOyEO E ROM DOUBLE PRECISION STATEMENT-, 
APPEARING- IN OTHER ROUT I Nt S USED .IN CONJUNCTION WITH THIS 

^ KS'JoSblE PRECISION VERSION MAY B^IMEO BY CHANG I MO THE 
CONSTANT IN STATEMENT 78 TO b.00-12 A^JN STATEMENT M 
1.5-10. THIS WILL PROVIDE HIGHER. PREUSION RESULTS, AT THF 
COST OF- EXECUTION tl ME •« t 



MK 1 

MK 1 1 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MKl 

MK I 

MKl 

MKl 

MKl 

MKl 

MKl 

Mfcl 

MKl 

MKl 



8300 
8310 
8320 
3330 
8340 
8350 
8360 
8370 
8380 
8390 
8400 
8410 
,8420 
'8430 
8440 
8450 
18460 
8470 
8480 
8490 
8500 



\ 501 
| 502 
! 503 
504 
505 



506 
507 



508 
.509 



510 
511 
512 
513 
514 
515 
516 
517 
518 
519 



520 
521 



522 
523 



524 
525 



10 

C 
C 
C 

20 

30 

C 

C 

C 

40 

c 
c. 
c 

50 

60 
70 



o0 
C 
C 
C 

90 
C 

c 
c 

LOO 

C 

C 

c 




IF1T=0 
N=M 

I6R?0 xj; ^ . 

IF (XCOMNH) i?10t40tlJ 

IF (MV 20,20/60 

v. 

StT error cere Tg 1 
IEK=i 

RETURN 

4tT EKKCk C( UF TO 4 

IEJR=4 r - 
GO TO -30 

SET cRRO* C^nfc TU 2 
IEP*2 

GO TO 30 * 
IF (N-36) 70,70,5Q 
NX = U 

NXX=N+I " 

N2 = l* 

KJl=N+l 

00 L = UKJi 

KT=KJi-hM 

CUF(MT) =XCtlF(L) 

S£T I'll \ \ A I. VAl IHS 

XO=.00!>00iOl « 
YUO. OiOJOlOi , * 



60 - 



ZtPC INITIAL V ALUt C'UNICh 
X = X 1 ' 

1,WEME.\'T'MTIU VALUES AND COUNTER 

X0=-LC.0*YC. ^ - 

Y0=-lO. J*X 

J X A NO Y TO LUtvUrfcT VALUE 

X = X'1 

Y = YJ 
IN«IN+l, 
GO JO 12 3 
I F IT= L 
XPP^X' * 

YPrl =Y 

EVALUATE PUL YN'j^IAL, ANO OFPIVATIVES 

ICT = 3 r, 

UX-3.0 

UYO.C 

V = 0.0 



MKi 
MKl 
MKi 
MKi 
MKI 
MKi 
MKI 

Mki 

MKI 
MKI 
MKi 
MKI 
MKi 
MKI 
MKi 
MKI 
MKI 
MKi 
MKI 
MKi 
MKI 
MKI 
MKi 
MKI 
MKI 
MKI 
MKI 
MK I 
MKL 
MKI 
MKI 
MKl 
MKI 
MKI 
MKI 
MKI, 
MKI 
MKI 
MKI 
MKl 
MKI 
MKl 
MKl 
MKi 
MKl 
MKl 
MKl 
MKi 
MKl 
MKi 
MKl 
MKl 
MKl 
MKi 
MKl 
MKL 
MKi 
MKl 
MKl 
MKl 



8510 
8520 v 
8530 
8540 
8550 
8560 
8570 
8580 
8590' 
8600 
8610 
8620 
8630 
8640 
8650 
8660 
8670 
&680 
8690 
8700 
8710 
8720 
<8730 
8740 
8750 
8760 
8770 
8780 
8790 
8800 " 
8410 
8820 
883-0 
8840 
8350 
8860 
8870 
8880 
8.390 
8900 
8910 
8920 
8930 
8940 
8950 
8960 
8970 
3980 
8990 
9000 
9010 
9020 
9030 
9040 
9050 
.9060, 
9070 
9030 
9090 ' 
9100 



65 



140 



560 
561 
562. 
563 



150 



160 



C 
C 
C 

170, 

180 
190 
C 



- 61 - 

YTO.C* 
XT*i.O 

IF (U)" 140,27** UC 
On 150 I-liN 
L-N-I+i 
T£.MP=C0F(U) v 
XT2=X*XT-Y*YT 

YT2-X*YT + Y*XT * 

o=u*Te.«p,*.xT2 

VsV + Tfi*iP>YT2 
k Fl = ! 

ux^x+f i**i *rf-MP ■ v 

UY = »JY-F I *Yr*TLrlP". 
XT =XT 2 * 
YT*YT2 
*Stm.} = UX*UX+UY* UY 
IF (SUMSQ) io0,^30, 16C 
DX= ( V*UY-U*UX ) / SUM^O 
X=XOX 

0Y=-(U*OY*V*UX J /SUMSO 
Y = Y*t)Y • 

IF I Ar*S(OY)+A35(»)X)-WOl-5) 210,170,170 

STEP lT6KATlt:N CHUNK* C 
1CT=ICT*1 

IF IICT-500) UOflrtCidC 
IF ( IFTT) 2U,190,riO . 
If ( I N- 5 J ,100, ? )J, 200 

StT £ Rkfjfc C T-i TP i t 



'4 ' 

•4*. 



564 


200 


1CP^=3 




36b 




GO TU 30 




566 


210 


00 III L=1,NXX 




567 




MT-\J1-L+1 




568 




T£MP=XCUF(MT) 




569 




XCOF I MT ) =C0F ( I J 




570 


220 


COF,(L) = TFMP 








I T FMP*N' 




572 




N=NX 




573 




- NX* I TEMP 




574 




IF f IF I T 1 250, 110, 


250 . 


575 


230 


IF (IFIT) 240,100, 


240 " f 


576 , 


240 


X-X»R ' 




577 




Y-YPR 




578 


2 50 


i f i r*o 




•579 




- IF (AeS(Y)-UO--'i< 


.% tv S ( X > J ?Hpf26Q. f 260 


5§0<i 


26Q 


AlrPHA=X*X 




581 




' SUMSQ=X*X+Y*Y 




5B2 




N=N-2 r 




583 




'GC TO 290 / 




584 


2 70 


,X=C.0 




585 




NX=NX-I 




586 




NXX=MXX-l 




587 


230 


y*o.o . 




588 




SUMSO=0.0 




589 




ALPKA=X 




590 




N=N-1 





MK1 9110 
MK1 9120 
MKi 9130 
MK1 9140 
MM 9150 
MKi 9160 

• MKI 9170' 
KKl 9180 
MKi 91^90 « 
*\l 9200' 
MM 9210 

-MKI 9220 
MKI 9230 
.^Ki 9240" 
MKI 9250 
kKl 9260 
MKi 9270 
MKI 4280 ' 
MKi 9290 j 
MKI 9300 
MKi 9310. ^. 
«K1 9^20 
MK 1, 9330 
MKI 9340 
MKI 9350 
MKI 9360 
MKI 9370 
MKI 9380 
MKi ^390 
MKi 9400 
v Ki 94l0- 
MKi 9420 > 
MKI ,9430 
MK j^9440 
MKI 9450 
MKi 9460, 4 
MKI 9470 
MKI 9480 
MKI 949fr 
MK JL 9500*- 
MKI <>510 
MKI 9520 
MKI 9530 
MKI 9 540 
MKI 9550 
MKi 9566 
MKI 9570 
MKI .9530 
MKI 9590 
MKI 9600 
MKI 9610 
MKI 9620 
MKi 9630 
MKi 9640 
MKI O650 
MKI 9660 . 
MKI 96J0 

S MK1 9680 
MKI 9690 
MKf 9700 



66 



591 
592 
593 
594 
595 
596 
597 
598 
599 
600 
601 
602 
603 



290 



• 300 



32Q 



3 30 

C 

c 

(y 

c 
c 
c 
c 
c 
c 
c 
c 
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• • 
CPF (*2 ) = COF ( 2 ) *ALPHA*C0F ( I ) 
IF (N.tO^O) GO TO 310 

' i/o 3oc 1*2; n 

COF<L+l )-COF( L*i KAlPHA*COHt)-SUMSQ»COF(L-ll 

*O0THN2)*Y 

R0pTR(N2)=X 

U2=N£>1 % 

IF (SUHSO) 320, 330, 3120 

s-uyso-o.o 

GO TO 3 10 ^ 

IF (N> 30,30,93 . * 

*ND 



SUUROUT 1*41 THLOT C l S b I NGLfe PRFCISION PLOTTING 5 DIFF 

VARIABLES » ASSU^l^G CONSTANT TIM£ INCREMENT BETWEEN TWCT 
EVENTS . 



SUPPORTING KLU-T I Ni NCNF 



MKl 
MKl 
MKi 
MKl 
MKl 
MKl 
v MKl. 

rtK'i 

MKl 
MKl 
MKl 
MKl 
. MKl 
MKl 
MKl 
MKl 
:RENTMK1 
>, MKl 
MKi 
MKl, 
MKl 
MKl 
MKl 
MKl 



9710 
97?0 
9730 
9 740 
9750 
,9760 
9770 
9780 
9790 
9800 
9610 
9820 
9833 
9840 
9850 
9860 
9870 
9880 
9890 
9900 
9910 
9920 
9930 
9940 



604 
60b 
606 
607 
608 
609 
610 
611 
6l2 
613 
614 
615' 
616 
617 
618 
619 
620 
621 
622 . 
62* 
624 
625 
. 626 
62 7 
628 
629 
630* 
631 
6*2 
6 3'3 
634 

636 
637 



10 



SUtf^JUTINf TPLjT H,l •M2iM6tH9,M0 l JXJ 
MPLICIT *CAL*4< A-h.^-Z) 

OIMFMbinvj ^8(JX), MO(JX), "C(JX), Ml(JX) f M2(JX) 
0 i M&.MS I >N LINFU1), VMJM(«*) 

iMFGF« PL,MI i^T,«LiSL,S^tbw,S6,SifS2 , 
READ [j f *C) PlffMf ST.HLfSl iS9,$0tSii$2tS6 
MXY^O, j * , 

MINO=0. Ji 
^IN2=0.0 
MlN8=0.0 
M1N9=0.0 
PHI3=0.0 
pMl2=0.0 
PHIfi=0.0 

PHI9=U.0 „ - , 

DO 10 I=liJX ' • 

(MlNO.GT.X^t I ) ) VINJ-MCtD 

(MIN2.GT.M2C U) *Ifc2*M2'( J ) . m " 9 

YlINa.GT.M8C I ) ) ^HvJj = M6( I ) 
CMIN9.GT .M9l I J ) >IN9=M9(J ) 
(A6b(M0( I) LCI. PH10) 
4ABSIM21 I) I.0T.PHI-2J 
( AbS ( MR ( I) J..GT.PHI8) 
(^i>S(M9t I I I.GT.PH19) 
CLNtlNUF 
JJ'JX . 

J JO=J Jjr'6 + 1 
JJ1=JJ*1 

kRI*€ (6,901 
WRITE (^,100) 
PHIO^PHI 0*Ah$( IN J f 
f*Hl2 = PHI2*A4S(^IN2> 
PHl3 = PHl8«-AB$m'f4S) 
PMI9»'PHI9*ABS(Mino) f * 



IF 
IF 
I F 
I F 
IF 
IF 
IF 
If^ 



PHI0= A^SCMJC I ) ) r 
PHI ?=£bSf*2( ! ) ) 
PHI 3= ABS( M8( I ) ) 
PHI9=A8S(M9( I) ) 



MKl 9950 
MKl 9960- 
MKL 9970 
MKi 9980 
MKl 9990 
MK110030 
MK110010 
MK110020 
MK110030 
MK110040 
MK110050 
MKl 10060 
MKlLOOja^ 
MK11O08O X 
MKU'0090 , 
MKilOlOO 
MK110110 
MKliai20 
MKI 10130 
MK110140 
MK1L0150 
MK110160 
MK110170 
MK110130 
1 MK110190 
MK110200 
MKi 102 10 
MK Iff 0220 . 
MK1L0230 
HK110240 ' 
MK110250 
MK110260 
MK110270 
110T280 



.7 



67. 



I? 1 



If* * 



Z7 



638 > 
639 
640 
641 

642 
643 % , 
644 
645 
£46 
€47 
64.8 , 
649 
650 
651 
652 
653 
65h 
,655 
6^56 
-657 
653 
659 
660 
661 
66^ / 

663 
.v«>4 

! 665 , 
^"666 
- 66J 
% 668 . 
.6.69 • 

670 
'671 

672 

675 

674 
"675 

676 
" 677 

678 

679 

680 ' 

681 

682 

664 
685 ' 
686 
687 
688 
689 • 
+ 690 

69 r 

~69£ 

693 

694 

695 

6«T6 
*697 



20 
30 



50 



-60 
70' 



dO 
)0 
100 
118 



>. - 63 - 

>0i i l»HO« I I*ABSCM1N0) 
M2(l)=K2(l)*AdS(MlN2) 
M4( l)=h8< I )+ABS<MIN8) 
NH(n = M9<! )t-ABS(M!N9) 



DO 2J I-ltJJ 
IF' M1N0.LT .O..o) 
JFMMlN2.lJ^0.O) 
IP (MIN8.U^3*0) 
IF (MIN9.U .0.0) 
M0( 1 )-MO(l )/PH10 
M2I 1 )=M2C I J/PHI2 
Mi I ) -M8I I l/PMId 
#491 1 )=M9(I JAPHl'/ 
UNO 3j; I = lf 9 
INtlMU ) = l 

I ff (6,110) CINUJMI)il*lf9) 
00 70 I = l\JJl 

if- (i .10. i ) on to i 

MXY=Ml(!-l) . * 

rPB=ie( i-i )*&o*i.u 

joo=M9( I-l)*6C*1.0 
IPG=&0( l-l >*oO*l.u 
IP2=K2( I~l)*6 0*I.C 
f 00 40 Il=lt56,5 
LlNc(Il)=8L 
DO 40 I 2=1,5 
I3=.l UI2 

If ( ll.fcQ.IPJ) LINt (ll)^SC 
(I3.EO.IP0) UU<r(|3)rS0 
(Il.CQ.IP2) ■UNCCIll*S2 
( I3.FQ.1P2) I < I3)=S2 
( I i.fcO. IPJ ) I INc ( H ) = S1 
( 13.F0.IP<i) tt tt' -t< 13 # >*«Sl 
( II.MJ.IP9) Ll'iH ll) = S9 
Uf:H |3I=S9 



IF 
IF 
"IF 
" IF 
•If 
IF 
IF 



(13. 



r 



c 0* I P9 J 
L<JM|NUF 3 
I INF( 61 )-PL 
|'IL = 1-1 
IF(l .eQ.I^O) Ll'. r M>*SO 
IF( I .E0.IP2) I IMP I t) = $2 
1F<1 .fOi IP8) LINK 1 J -S 1 
!F(1 .EU.IP/9) U'.F(1)=$9 , 
IF( lPC.FQ.f l)Ufc£(61)='iG 
IF(IP^.e0.6l)UI^fc(61)=b2' , 
IF( 1P8.Q0.61)LINHM)=S1 

IF(lP9.F0.6i)L!Nr (61)=S9, , .»...» n . 

I'f-dP^.NC.r.Q^.l^fi.Nt.l.nKWPO.NC . 1. OK. IP2.NS:. I )UNU U^U 
WR ITE( 6,120 )MXY,f LI NHKK),KK = 1,61 ) \ 
lFU.fO.JJl) GO TO 70 # \ 

GONTINUF , 4 " , ' " . ^ 

00 60* 1 1=1 06,5 - - t , * , * ' 

DO 60 12=1 • 5 / fc 1 ' - , » 

13=1 1*1 2 

LlNF(13)=RL_ - . 

C0NT1NUF . 

CONTINUF ' " ^ 

WKITE(6,130)( lMJM(l)..I = 1.9l ■ A . 

. WRlTF(6,lt>C)Pt,M.St VBL,SL,S9,S0,S1,S?,S6 ^ ' ^ 

WRIT£(6,160) 7 1 " . 

WRITF (6,140) 1 
FC'PMAH UAf) * ' 

FORMAT ( 'PI y ^ m ' 

FORMAIOSXt'P.fcLATlVt f)£r.$lTY» ) 

FOF MAT (i7X,9(2X.»^SU» < 1 »/il4XtlOC«* -•)••♦•> Q 



MK110290 
MK110300 
MKU0310 
MK110320 
MK110330 
MK110340 ' 
VK110350 
MKl 10363 
MKU0370 
MK1103S0 
MK 110390 
MKl I 04.00 
MKl 10419 
MKU0.420 
MM 10430 
MK I 1 0440 
MK110450 
rtKU0460 *" 
MK11047C 
MK110480 
MKl 10490 
MK110590 
MKl 10510 
MKU0520 
MK110530 
MKl 10540 ' 
MK110&50 ' 
MKl 10560 
HK110570 
MK1105^0 
MK110590 
MK*10600 
^K110610 
MKU0&20 - 
MKl 10630 
MK110640 
MK110655 
MKl 10660 
MKil0670 
MKi 10680 
MK1106?90 
^Kl 10700 
MK 110710 
MK110720 
^ MK113730. 
MK1L0 740 \ 

MKU0750 ; 

' MK1107<>0 " 
Nj^l 10770 * 
MKl 1 07 80 
110797) 
MK 11 0930 

MKIIQ320 — 
M^l l^^30* * 
MK110840 

MKiioaso" 

MKl 10^60 ■ 
MK110870 " 
KK11O880 



64 



■ 69d 
699 
700 

'701 
702 



/03 
704 



120 
130 

150 



FORMAT (IX, •TJME , ,lX.h7.3, IX 9 6 LAI) 

F.OFMAT ( 1<*X, iOt 1 ♦ « ) i/t 16X,9(3X, »0. • ,11) I 

"*tJP MAT ( • 1 • ) . 
FORMAT ( iSx, • INPUT CHA U AC T I k • 1 1 1 A l f / ) 



MK110890 
MK110900 
MK110910 
MK 11 0920 

16D^ F-CMATP S/.lbX, 1 P - P0W£K • ,/ , 1 , « L -MK11D930 

< L vG. OF P0WE*» ,/f 15X, ■ ij - P*FCU«$QR DfcNSITY* ,/,15X,» R - REACT1VMK110940 
. <ITY» 9 ///t U>X,'NEbAMVt» VALUE PLUTTEO WITH THE AXIS IN THF C ENTER 1 t MK 1J. 0950 
* <////) * MK110960 

STCP * MKU0970 

£^D „ * MK11093^ 



//DATA 



f 



\ 
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****************************** KINET ICS ""M&?)uLt: 1 ******************** 



- **;*************«*** 

INPUT PARAMETER 
******************* 



\ 



* TYPE OF FUEL 

NUMBER t)F O6LAYF0 ,N6utRON = 
NEUTRON GENERATION . T !ME = 
TYPE OF REACTIVITY * 2 



O.OOOiO '(SEC) 



RO = * 



0.25>0 



32 = 



5.00000 



i.000 (SEC) 



TOTAL J1M6 USED = 
OUTPUT UPTIQN = 2 
INITIAL RELATIVE FLuV— 1.300 



TIME INTERVAL = O.OiO 



********* gNO OF IMPUT DATA ********** 



SIX GROUP LAMBDA AR£ 0.0126 0*0-301 0. 1240 ^' 0, 3250 1.1200 
SrX GROUP bFTA AGC 3.0'jfOV 0.00080 0.00057 0.00088 0.00023 0 



THE WER AGE BETA AND LAMl) A 

B< L)= <D.J00°0 LAMdDA( 1) 

Bt 2)= 0.001*5 LAMdDA ( 2) 

B( 3M ty00035 s LAMBDA ( 3) 

COEFFICIENTS OF INHOUR FORMULA 



0,02636 
0.19854 
1.39571 







i 


• » 


M 


11 


= 


-O.OOOOE 00 










Al 


2) 


3 


0.30516-03 


A( 




= 


0.36066-02 


A< 


4) 




O.Z8<>26-02« 


A( 


5) 


= 


0-l?)06*O3 
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THE LAfcGfcST t lubN-VALOt 



NO TIME • R ^POhcR ' 3 OPOUP D6LAYE0 NEUTPGN 

< ) (S6C>\ <*) ( ) (RELATIVE) 



1 


0. 300 


0 f 00— 


1/ . J J r 


01 


0.3411r 


, u3 


0.7303F 


C2 


0.25 156 


ci 




T HE* G-MATK I X 












* 






•o. 


7634E 00 


0;2ilOE 


— 0 3 0 . 


174 




.12236-01 








0.899 1E-01 


0.999 7fc 


no > • 


:oc 


CE 00 J 


.00006 00 








0. 


1448E 00 


0. Ov ; r 'F 


J J 'J. 


we 




.OCOOf OQ 




- 




' 0.34866-01 


o .00 JCE 


00 0. 


OOCuF 00 0.O661F 30 








2 


O.^QO 


0 . } 


0 . 1 0 Ovj h 


01 


0.341U 


03 


0.73C36 


02 


0.25156 


01 


3 


0.010* 


0.01 


iy . i v y j l 


CI 


.341H 


03 


0.73036 


02 


. 0.251'^E 


01 


4 


0.020 


0.02 


0. lOOdC 


'H 


0.34116 


03 


0.7303E 


02 


0.2515f 


01 


5 


0.030 


j. 04 - 


O.lCl^t 


\j i 


0.341U 


05 


0.7303E 


0? 


0.2515C 


01 


6 


0.040 


o.os 


0.1024^, 


01 


0.3411f 




9.73036 


u2 


0.25166 


01 


7 


0.U50 


0. 06 


0.1u33Q 


01 


J.3411L 


03 


0.73C3E 


02 


0.25176 


01 


8 


0.^60 


u . 07 


0.1043*" 


01 


0.341if 


03 


0.730^6 


02 


0.25186 


01 


9 


0.070 


0 . J r > 


o;io5^r. 


c.i 


0 . 3'4 1 1 F 


C3 


0.7304t 


02 


0.25196 


01 


10 


0.330 


0.1 ) 


0 • I066t 


CI 


0.3411L 


C3* 


0.73056 


0£ 


0.2521T 


01 


11 


0.090 


0.11 


Q.1078E 


01 


6.3412F 


03 


0.73066 


02 


0.2523F 


01 


12 


0 . 1 0.0 


* 0,12 


o. IT9U 


01 


0.3412*- 


03 


0.7307E 


02 


0.25266 


01 


13 


0.110 


0.13 


0.11036 


CI 


0.3412E 


03 


0,730<*E 


02 


0.25296 


01 


14 


0*120 


0. 14 


i>.Tll6£. 


01 


0.3412F 


03 


0.73106 


02. 


0.2532E 


01 


'15 


0.130 


CIS 


0.1129L 


tl" 


0.34126 


03 


"0.7312E 


02 


0.25366 


01 


16 


0.140 


0.16 


0.11431 




0.34121- 


03 


, 3.7314E 


02 


0/2540E 


01 


tl 


0.150 


0. 17 


9.H5o£ 


CI 


0.3412E 


c3 


X). 73166 


02 


CV£5456 


01 


18 


0.160 


0. 13 


11696 




C.3412F 


03 


0.7U86 


C2 


0.25506 


01 


19 


0.170 


0.4^ 


0.1l82h 


01 


0'?341 2E 


' 03 


0.73206 


02 


0.25^56 


01 


20 


0. 180 


0."2J 


0*11956 


CI 


• 0.34136 


C3 


0.73236 


02 


0.25616 


01 


21 


0.190 


0.20 


0.l20fL 


CI 


0.^413F 


03 


0.73266 


02 


0.25676 


01 


22 


0.200 


0.2:1 


0.12206 


01 


C.3413E 


03 


0.73296 


02 


\ 0.25746 


01 


23 


0.210 


0.22 


0. 12326 


01 


0.3413* 


03 


0.73326 


02 


0.25816 


01 


24 


0.220 


0.22 - 


0.12446 


ol 


0.34136 


03 


0.73356 


02 


0.25886 


01 


25 


0.230 


0.23 


0.12556 


01 


0.3414E 


03 


0.73386 


02 


0.25956 


01 


26 


0.240 


0.23 


0.1266F 


C 1 


0.34146 


03 


0.73426 


02 


0.26036 


01 


27 


0.250 


0.24 


0.12766 


01 


0.3414F 


03 


*0.73466 


02 


0.26116 


01 



V 



67 - 



i ™ 


0.260 


{ 29 


0.27a 


« 30 • 


0.280 


! 31 


0,290 


32 


0.300 


33 


0. 310 


34 


0.320 


35 


0.330 


36 


0, ><»C 


37 


0.350 


36 . 


0.360 


39 


0.3 70 


40 % 


0.380 


41 


<0.39U 


42 


0.400 


43 


0.410 


44 


0.420 


45 


0.430 


46 


0.440 


47 


0\4 50 


48 


0.460 


49 


0.470 


50 


0.<t80 


51 


0.490 


« 52 * 


0.500 


53 


0.510 


54 


0.520 


55 • 


0.5 30 


56 


0.54o 


57 


0.55U 


58 


0.5oO 


59 


0.37t 


60 


0.390 


61 


0.59O 


62 


0 .oOO 


63 


0.610 


64 


0.620 


65 


0.&30 


66 


0.640 


67 


O.o50 


68 


0.66/0 


69 1 


0.670 


70 


0.680 


7 L 


0,o90 


72 


0.700 


73 


0.710 


74 


0.7?0 


75 


0.730 


76 


0.740 


77 


0..7 50 


78 


0.760 


79 


0.770 


80 


0.780 


81 


0.790 


82 


0.300 


83 


0.310 


84 


0.320 


♦ 85 


0*630 


\ 86 


0.340 


1 87 


0.850 







J. 24 


0, 


123ol 


CI 


0.3414c 


03 


0.7350F 


U2 


0.2619F 


01 


0.24 


0 .1295C 


Ci 


0.3414E 


03 


0. 7354F 


02 


0.2628b;* 01 


0.25 
*0.25 


0. 


1303b. 


oi 


0.3415F 


0 3 


0.735 7 8t 


02 , 


0.2637-E 


01 


0. 


13 11L 


* 1 


0.3415T 


03 


0. 73626 


02 


0.2646? 


01 


0.25 


0 • 


13186 


01 


0. 3415b 


03 


0. 7367E 


02 


a. 2655E 


01 


D.2i 


0. 


1324b 


Ll 


0.3416E 


0 3 


' 0. 7371E 


02 


0.2664f 


Oi 


0.25 


0. 


1329F. 


01, 


0.3416b 


03/ 


-0. 7376E 


02 


t). 26 73 b 


01 


0.25 


0. 


1333E 


CI 


0.34l6h 


03 


0. 7380F 


02 


0.2682b* 


01 


0.25 


0. 


1336b' 


C 1 


0.34i6r 


0 3 


0.7385E 


02 


0.269H 


01 


0.25 


0 . 


13381 


oi 


0.34 1 7F 


03 


0.73°0F 


02 


0.2701T 


01 


0 .24 


0. 


I339f 


01 


0 .3*1 76 


03 


0 • 7394c 


02 


0.2710F 


C I 


0.24 


0. 


13 39E 


f'l 


0.34176 


03 


0. 73996 


02 


0.27 19F 


01 


0.24* 


c . 


1338E 


ui 


0 . 3 4 1 8 F 


03 


0. 7404E 


C2 


0.2728C 


01 


J. 23 • 


0. 


1336E 


Ci 


0.3413f 


0 3 


'0. 7409E 


02 


0.2737F 


01 


0.23 


0. 


i: j 33E 


01 


0.34l8t 


03 


0.7413B 


02 


0.2746F 


01 


0.22 


0. 


1329b 


Oi 


0.3419C 


03 


0.74186 


02 


tf.£754F 


01 


0.22 


0. 


1324b 


ci 


0.34J.9F 


03 


0.7422E 


02 


t).2762F 


01 


0.2L 


o * 


mat 


Ci 


0. '34 191 


03 


0.7427E 


02 


0.2770F. 


01 


0.20 


0 . 


I3li c 


01 * 


0.3419F 


03 


0.7431{: 


02 


- 0.2778B 


01 


0.19 


0. 


I j03l 


CI 


0.342C? 


03 


0.7435E 


02 


0.2785b 


01 


J. 19 


0 . 


I295r 


CI 


0.3420C 


03 


0. 7439F 


02 


0.2792E 


01 


0. 18 


0. 


1205b 


01 


0.3420C 


03 


0. 7443E 


02 


' 0.2798F 


0 1 


J. 17 


0. 


12 75b 


CI 


o v 342oe 


03 


0* 744 7E 


02 


-^.28046 


0 1 


u .16 


G . 


1264b 


01 


0.*3421F 


03 


t 0. 745 Ih 


C2 


OV2010F 


01 


0.15 


•0. 


l'2 53t 


Ci 


* 0 . 3 4 2 1 C 


C3 ' 


0 . 7455F 


02 


0.28 1.5E 


01 


0.14 


0. 


1240b 


C 1 


C. 3421F 


C3 


0. 745 8 E 


0? 


0 .2819b 


01 


,0.13 


0. 


122-3*: 


'Jl . 


[} . 34 2 1 b 


03' 


0. 746 IE 


02 


0.2324F 


01 


0.12 




I21r>b 


Ll 


0.54$2b 


03 


0 • 7464-b 


02 


0.28 27b 


0 I 


-w. 1 I ' 


0 • 


120lt 


Ci 


0.3422E 


03 


0. 7467b 


02 


0#2830b 


01 


0. 10 


• 0. 1 I3db 


0] 


0.3422b 


03 


0. 7469b 


02 


0.2833b 


01 


j . : 3 


C . 




01 


0.34221 


ul 


0.7472T 


02 


0.2835f 


01 


D .c r 


0. 


1 16 


01 


0.3422b 


C3 


0.7474b 


02 


0.2837F 


01 


0 * Oo 


0. 


i I45h 


01 


0. J422fc 


03 


0.7476F 


02 


0,.28 3<8f 


01 


0 .65 


0. 


11311 




0.3422E 


C3 


0.7478? 


02 


0.2838b 


01 


O.t't 


c ♦ 


nut- 


01 


0.3423b 


03 


0.-74796 


02 


A.2Q38E 


01 


J.j2 


v. 1 • 


l iC2t 


Jl 


0.3*23c 


03 


0.74*6 IE 


02 


0.2838F 


01 


U. r 1 




ICSof 


Ci 


0.3423b 


C3 


* 0.7482b 02 


0.28^7b 


01 


-a. oo 


0. 


1074C 


> 1 


0.3423E 


03 


, 0.7463b 


u2 


0. 2836E 


01 


-0.C1 


0 # 


106 3b 


"r i 


0."3423t 


03 


0.7463F 


02 


0.2834F 


01 


-0.C3 


c. 


1046b 


Ci . 


„0. 3423b 


03 


0.7484C 


02 


0.2832b 


01 


-0.04 


0. 


1034F 


01 


0.3423b 


03 


0. 7464E 


02 


0.2829b 


01 


-0.05 


0. 


iC2Ufc 


CI 


G.343 3F 


03 


0.7484E 


02 • 


0.2826E 


Oi 


-0-.06 


0. 


IC J7fc 


01 


0.3423fc 


03 


0.7434E 


02 


0.2822b 


01 


~ J • C8 


0. 


9947b 


CO 


0. 3423F 


C3 


.0.7484F 


02 


C.2818E 


0 i 


-0 .09 . 


,G.922oE 


00 


0.3423E 


03 


0. 7483b 


02 


0.28 14b 


01 


-0. 10 


0.9703E 


c:, 


0. 3423b 


C3 


0.7483b 


02 


0.2809b 


01 


-0.il 


C. 


95956 


C j 


0.3423b 


03 


0.7482E 


02 


0.2804b 


01 


-0.12 


0. 


948?ir 


C j 


0.3423E 


C3 


0.7481E 


02 


0.2798b 


01 


-0. 13 


0. 


9330!: 


ca 


0.3423C 


03 


0.7480E 


02 


0.2793£ 


01 


-0 .1*4 


0.9279E 


CO 


• 0. 3423E 


03 


0.7479F 


02 


0.2787b 


01 


-0. 15 


0.9182b 


Ou 


0.3423E 


03 


0.74776 


02 


0.27801 


01 


-*0 . 1 f> 


J • 


9089b 


00 


0. 3423E 


0> 


0.7476E 


02 
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REACTOR DYNAMICS MODULE, RD-2 

... to 

REACTOR KINETICS WITH FEEDBACK 

f 

2.1 Object of Module 



The object of this module is t^ . , , v 

(1) Examine the temperature feedback mechanism of a PWR ai*d 

(2) Solve the one delayed neuron model with' temperature feedback for 

. a step insertion and a ramp insertion of reactivity. « 

The time dependence of a reactor, taking the feedback mechanisms into 

account, is relatively difficult. We will consider a PWR core with a two path 

feedback. The reactivity is diminished as i^the temperature of the fuel increases 

♦ 

due to the Dopf>ler broadening of the resonances . 4 This feedback is instantaneous 
since the temperature increase follows the power generated immediately. The 
secqnd feedback path is that of the moderator temperature coefficient. As the 
moderator temperature increases, the number density decreases and the neutron 
mean free path increases so that "leakage increases and reactivity decreases. 

WTe i^ll be concerned about the stability of the reactor; to a limited degree, 

-The dynamic response depends upon the magnitude ,of the temperature coefficients 
as weU|as that of the signs. For a given reactor design* i.e., a given life- 

'time ^Bkpow^r level, the reactor «&y or may not be stable for a given set of 
reactivity coefficients. 

The core region is the only one of interest in this module. The r-pst of 
the primary loop as well as the secondary loop is treated in an overall dynamics 
module for a PWR. " „ 
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Also, all of our analysis will be f und amen tal' mode analysis. The physical 

i 

phenomena are taking place so slowly that the higher "harmonics of the flux 
distribution are all dying out so rapidly that we only need to consider the 
lowest'or fundamental mode. 

The thermal analysis really should proceed by the solution of the space- 
titoe heat <$cmduction equation. This is a very complicated procedure- and would 
also mefcn that spatial effects of the kinetics equations should be taken into 
account ^ We, therefore, will assume only a lumped parameter model and will 
obtain the time dependence of a reactor which is r^lly one with the average 
properties of the reactor under consideration. 

The program name is FUMOTEM which is an acrorfym for 1 "Fundamental Mode 

Kinetics with Temp erature Feedback." 

' v . & 
There are four types of reactivity inputs that the program can accommo- 

'j 

date with NRO = 1, 2, 3 or, 4 respectively,: 

1) P Q (t) = P Q o < t < t r • * 
k • = o otherwise 

2) p o (t) - p o (l + a t) o^t<t r 

- * p^(l + a t r ) otherwise 

3) P„(t) 5=1 'Qt cos a t 



o 



or 



4) P n (t) « p sin, a t ' 

The feedback reactivity is taken to be rerp for times t^ f One *c§kr feed 

in a t* greater than the calculation time however. So for*t' > t , the only 
r . • r • , 

reactivity is due to feedback effects. , 

• '''#■ 



2.2 The Feedback Model 



We will use the one group delayed neutron model to describe*the core 
neutronics. The kinetics equations are 



fiX£U £*£1jL! P(t) + x Q (t> 



(2.2.1) 



and 



^l-.f P(t) -XQXt) 



(2.2.2) 



where 



.t 



P(t) = The total reactor power (Megawatts) 
Q(t) '« Thet power equivalence of the delayed 
neutron precursors (Megawatts) . 



Now we let AT be the deviation of the spatially averaged moderator 
M ^ ^"^ 

temperature from its equilibrium value, i.e., 



AT M (t) - T M (t) - T MQ 




(2.2.3) 



and similarly for the fuel temperature T we have 



AT F (t> = T F (t) - T FQ 



A2.2A) 



wh^re\r and T * are the equilibrium moderator and fuel temperatures re- 

, ^MO- - FO •* 

spectiveiy. Also, we let and be the moderator and fuel temperature * 
coefficients of reactivity. Then (1) ^ ' 

P(t) - P 0 (t) + OjjAT M (t)^ F A F (t) . 9 (2.2.5) 

"5 •• 



I. 0 

jERjC 
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where generally; and a p will be negative or at least their sum is negative. 
The temperature coefficient of reactivity for the ^moderator is 



3d. 



1* 9k eff 



and 



"*M " 3T " k 3T M 
M ef f M 



= 1 9k eff 

aF=3T F 0 " k eff 3T F 



(2.2.6) 



(2.2.7X 



The thermal analysis of the core must now be considered "and Connected to 
the neutronics. The heat generated depends upon the fission rate or the. power. 
We will look at the fuel temperature averaged over the core, "as well as an 
averaged coolant temperature. We will ignore the cladding of the fuel pins. 

The heat balance equation for reactor fuel is * 



Rate at which the 
internal energy of 
the fuel changes 



Production rate 
of the energy 
in the fuel 



Rate at which 
heat is conducted 
out of the fuel 



or 



tfc 

where 



dX p (t) 



PF C F V F -It" = P(t) " 4,Tk F L F (T o " V N 



(2.2.8) . 



P F ' - density of 'fuel (lb/ft 3 ) 

C_* = specific heat of fuel (Btu/lb- 6 F) 
r 

V,, = volume of fuel* (ft 3 ) 
F 

* 

T« = average temperature of the fuel (°P) 

P(t) = total jJower of reactor (Btu/hr) 



= thermal conductivity of fuel (- — ^ tU OT J 
F . hr-f t- F 

Lp = length of fuel pin (ft) 

T - centerline temperature of" the fuel (°F) 

o 4 v 1 

- • f y' 

T temperature of the fuel pellets at the outer edge 

£ (pellet - water interface) (°F) 

N = total number of fuel pins in the reactor. 

i 

The expression for Equation (2.2.8) was obtained from El-Wakil, ,! Nuclear Heat 
Transport", page 123, equation 5-48. 

The^ temperatures ar6 averaged over the fuel pins. If we number each of 
the pins in the core, then the centerline temperature is 

• , T o = i ( V^ T o2 + •■• • + V 

where T is, the centerline temperature for the ith pin. T is defined 



^ 



similarly so that it As thp spatially averaged temperature of the "average 
fuei pin." \ s> ' s - 

Equation (2.2.8) contains the centerline temperature and the temperature 

at the edge of the pellet T vhich -must be eliminated. We assume that a 

R 

parabpli* temperature distribution holds even in the transient situation (really 
the transient heat conduction equatiorf holds here) so that for one fuel pin we 



have 



-2 

T(r) - T a . - ' 
F 



where 'T(r) is , the temperature of the fuel pin a distance r from the, c enter o f, 

K * \ J 

the pin. the average fuel temperature* is then 



•x 

N 9 



1 



T (t) « / rdr 



T - 



P(t)r" 
4k F V F 



= T 



P(t) ' s 2 
o 8k F V F *p ' 



(2.2.10) 



Eliminating the centerline temperatures, we have, using (2.2.10), 



t p (t) = T R (t) + 



R^ 2 P(t) 
~ 1 W 



.(2.2.11) 



Equation (2.2.8) can now be written as 



dT F (t) 



p F c F v F -ir^ = p(t) - 47T W 



p(t) *v 

T (t) + T (t) 



. =|l/2jp(t) - 4itkpNLpTp(t) + 4itkpLpT^(t) * (2.2,12) 



The wall temperature of the peLUtfi^T^t) is connected to the coolant; 



temperature * since 



(2.2;i3) 



where f A is the total,. area of the fuel and h is the heat transfer coefficient 
F 

2 

for the fuel water interface (Btu/(Ft ,hr-°F)). 

The energy balance for the water in the core is 



[Heat stored in 
water in core 1 



Heat conducted 
in f rott fuel 



Heat transferred 
out of reactor core 
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or mathematically, 

dT M (t) . 

"mVm Hit- = *• V* V«rW - V t)J + Wmi - S^V -< 2 - 2 - 14 > 

where m is the mass flow rate of the water through the core and T (t) is the 

M M 



average moderator temperature. Also, assume C - C • 

Ml MZ« «s 

' We will assume is an- input, as is T , the moderator inlet temperature. 

Thp moderator outlet temperature is related to the inlet and average temperature 



as 



T M2 (t) = 2T M ( ° " T M1 
so that Equation (2.2.14) becomes (using Equation (2.2.13)) 

' c v dT ^ (t) = 2irtLL- • + 2^0, £T -«T (t)] 

' * M M M dt TFT \\ M Ml M 

. • - ■ ... • v 

~> . P(t) + 2 - 2 ^„(t)^ (2.2.15) 



This can be rewritten as 
dT 



t . JX V» 

2m 



T \ * 

2 

The area of the fuel i^ 2NL ttR and the fuel Volume is ttIL, L N. 

F F x V ¥ 

In summary,* the equations we must solve are Equations (2.2.1), (2.2.2), 

• / 

(2.2.12) and (2.2.16), the last two of which are written J . / 

dT_ . ^. . . 

• ' p F c F v F^r = -(1/2)p(t) " ^WV 0 + HWV^' (2 - 2 - 17) 



d V c > . jssl_ _ ^ T (t)+ _2_i i • 



dt 



put these equations, into matrix form we have 



(2.2.18)- 



P(t) 



Q(t) 



T F (t) 



T M (t) 



P-B 
A 



2 >F C F V F * 



P C V 
M M M 



P(t)~ 



Q(t) 



T F (t) 



.0 



4irKpLpN 

pCV„ 
F F F 



' 0- 



0 



4itK f L f N 

Vm 



M M 



(2.2.19) 



'9 



where the h is dependent upon the temperatures themselves. Thus, Equation 
(2. 2. .20) is non-linear. 

The linearization of Equation" (2.2.19) (or (2.2.20)) can be achieved 
assuming that we can look at changes about some operating point. Expanding 
about the operating point we have 

P( t ) = P° + AP (2.2,21) 

Q(t) = Q° + AQ (2.2.22) 

T F (t) = TJ + AT p <> (2.2.23) 

T (t) = ' T° + AT (2.2.24) 
M M M \ 

We assume that P* etc. are independent of time so 

* * 0 

; -d A* ' . J 

r — - — = (A° + AA) (<t>° + A<f>) + B 1 , 

^ dt = = r- - 

or neglecting the AA • A<j> term we have 

— A6 = A° 4>° + A° A<{> ,+ AA ^ + B . 
dt 1 =* 1 ■ - = ~ 



or 



— A6 = A° A<f> + AA6° 

dt 1 s - -* \ (2.2.25) 



since 



A° £° + B » 0. 



• 




Ierjc : .. - r v;^-i/.-.'--v:- 




X 
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Writing Equatioif. (2.2.25) out explicitly, we have 



d_ 
dt 



AP 



AQ 



AT, 



AT. 



M 




or 



d_ 

dt- 



AP 



AQ 



AT T 



gERJC 



o 

A. 



-X 



>f c f v f 



4Trk F L F N 



>f c f v f 



-2 



. 0 • 0 



p M V M 



AP 



Aq 



AT, 




AT, 



M 



VW t m : * 



£. 
A 



0 0 0 



0 0 0' 



0 • .0*0 0 

I 



0 ' 0 0 



L MJ 



AP 



-X 0 



2 ?f c f v f 



Wf P F C F V F 



-2 



'mVm 



0 

.88. 



>M V M 



! * 



AQ 




(2.2.26) 



# • 
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Equation (2. 2. 26). is a linearized version of Equation (2.2.19) . This 
equation is soluble by the ordinary method of eigenvalues. In matrix form, 

9 „ 

this can be written as * 

dA 4 



dt 



A' A£, " (2.2.27) 



and A 1 is given in terms of the equilibrium values. 

The system of Equations (2.2.19) is highly non-linear because of the 
first element of the mattix A. The method discussed in Kinetics Module 1, 

developed by Hansen, is still applicable even to non-*lin&ar systems. 1 There is 

f 

a slight modification, that we must discuss and this will be outlined in a 
later section. 

•The overall heat transfer coefficient is also ffelatiVely difficult to 
obtain so we discuss this and related parameters in the next section. 

' ' * / ' ' - 

Problem 2.2.1 

Show. that if T is the temperature in the ceriter of a fuel pin and T is 

o . • * - t- n 

the moderator temperature,' then . 



T 1 T , P(t)a 2 ., P(t) 



1 ? b 1 
o *M • 4k_V F • 2V p J k clad a li^b 



The- inner radius of the clad. is "a" and the outer radius is "b.' 



I 



2.3 The Feedbaqfc-Parameters 




There are ^various nafanieters that must be obtained in order to numerically 



solve Equation (2.2.\9). We list these parameters in order: 



M 



The moderator (coolant) temperature coefficient of 
reactivity. 



2 ^£? - " 1116 tem P erature coefficient -of reactivity for the fuel. 

3. Cp F -^The total heat capacity of the fuel. ' 

4. k^, ^- The^ thermal conductivity of the fuel. 

5. h T - The overall convective'heat transfer coefficient of the 

fuel to numerator surface. 



The flow rate F, the fuel volume V , etc. are all parameters that can be 
obtained* from specifications of a particular teactor type and we shall not 
discuss these any further. 

Wg first obtain a relation for c^. The Final Safety Analysis Report of 
^nuclear power plants generally have a curve of ^ as a f unction^ of power level 
for a given boron insertion, and a critical rod insertion.- The curves generally 
follow the pattern that ct^ is constant from about 10% to 60%' power and then 
decreases rather dramatically between 60 and 100% power. a M is assumed to b^ 



negative. 



-4 Ak 
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.4 - 

.2 - 
.1 -- 



50 

Power % . 



H h— f 



100 



We will therefore take 0^ to be an inpufe-cpnstant but an' accurate analysis 
would necessitate a knowledge of » the variation of with power. 



The OpCP) is also rather difficult to calculate. ^To do this, we use 
Equation (2.2.^),- i.e. 

) 



_V 3k eff 
^ k eff 9T F 



Nw if we use the relation 



(2.3.1) 




k eff = fe L f L th (2 : 3 - 2) 



£n(k gff ) - £n(n fe L f L th ) + £n(p) 

and if we assume that the resonance escape probability p is the only factor 
which changes with the fuel temperature, then 



dk 




„ _ L_ eff 1 dp_ ft 3^) 



A standard expression for the resonance escape probability is (2) 

~ir v 



* SM-M ' . / o O AN 

p ^ e * (2.3.4) 

aJ^T the resonance, integral R(Ty) is given by the" empirical relation 

\ * - . 

R(T p ) - RG )[T» Y(^-^)l 4 ' (2.3.5) 

— * ♦ • 
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where is the Equilibrium temperature of the fuel. 



We 



now relate ff p ff to a in a way s'o that* the a can be calculated. Let 



then 



and 



£n p = - a 5 R( V , _l_ (2.3.7) 

) 

dR(T ) % R(U y 
* ' < dT F 5 2^1 ■ - 



IF we substitute T F0 into p we gdt 



or 



For an arbitrary value of temperature Ti, we obtain 



Dividing Equation (2.3.9) by Equation (2.3.10)- we. get 



An .u 1 , * a c R(T„). ■ ' v (2.3.9) 



£n = a 5 R(T p ) . ' (2.3.10) 



\ . ■ • 

Inserting Equation (2.3.11) into Equation (2.3.8) we obtain 



y ' Y3 5 • * n * (T F0 ) 



' • ... ?j* - 2 /r 15 v F0 ' J 2^: v / . * n p <V 

4 * F, » ^ 

Equation (2.3.10) yields ■ , - 



a^ R(T p ) % - £n p(T p ) 



so 



*» P(T F0 ) 



F ---t-i- My ^ 

F 



_ ' 1 (2.3.12) 

~ 2^ * n P (T F0 } .^ y 



In our modple, we will assume |hat y and p(Tp Q ) are read in. 

The^ third item of discussion -is Cp^, the total heat capacity of the fuel.' 



It is obvious that 



• . CpF " V F C ?Y * P F ,'(2 .3.13) 

where C* is the specific heat capacity of the, fuel, i.e. the units are 

Btu/lbm-°T% p-, is the density of fuel and V- is the volume of fuel. There is 
r v r 

a -small variation of C pF wi/th temperature but it is small for U0 2 in the 

Vegion of interest. The value that will be of interest for us is .0.0590 Btu/ib«-°F. 

\ * 

The thermal conductivity of the fuel is a function of the temperature. 
El Wakil (3). lists values of as a function df T p , We utae these values to 
obtain a polynomial regression df 'the' v k with ihe T F . 
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The last quantity- that we discuss is the overall convective he v a^ transfer 
* % * 

* Coefficient tip- The heat transfer coefficient is defined by Newton's law of* 

cooling* The relevant. relation is ~ 4 9 

where A^, is the fuel area. There are many factors which influence h such as 



* i) the temperature of the system 

ii) ' the heat flux > ' 

. ** > .' 7 

iii) the physical properties of the moderating material % 

* iv) the geometrical shape- of the* cooling surface « 

^v) the flow rate of tt& coolant 

Irv the PWR, the coolant flow is turbulertt. Therefore, to obtain the heat 

* • - ' • s. ■ - 

transfer coefficient* h^ we assume that • * 



" D ' " ' ' (2.3.15) 

e ■ ^ 



. \. #here % > * , • * 

*y * * • * ; * ' - . x 

* D s The equivalent diameter o£ flow channels through the-fuel 

' V v e v ^ * # . . ^ % \ 

• • rod b'undles ' \ ' * * - • 

Ic, =v The thermal conductivity' of ,tfre moderator f ^ " ) ^ J 

- * • ... , . 

C - The Colburn corfcection fact6r for .fluid flow paraMel to the 



tube bundles 
Re ** Reynolds number 



Pt * Prandt.1 number of the 'coolant in the core 
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Equation (2,3 / . 15) was a correlation recommended by Weisman. The Prandtl -number 



is taken from- the 1967 'ASMfe S.team Taliles, 



The < equivalent diameter of the fiow channel , is Ahown in Fig. 2.3.2. If 
^e let "a 11 be the area o^f^he flow channel and P be ttie Wetted perimeter, then 



w 



4A 

D e = F 
w 



(2.3.16) 



\ 



Foy a typical PWR, D, the diameter of the fuel rods is 0.03583 ft and the pitch 
(the distance between cun'terlines of fuel rods) is designated as P^ The 
Colbum correction factor is an empirical relation and is' * - . 

P 



V-c ="0.042 ~ -0.024. 



(2.3.17)' 



So we denote the- heat transfer coefficient, which varies with powef and 
therefore with time, as-h T <t)*and tmpequilibrium value as h(0). UsJDn'g 
Equatibn (2.3.15) we have — - 



ox 



***""' t y *•» * tyt) Re 7 f . 



!*T W y 0 ) ^ 0 - 8 <0) ?r l/3 (0) 



>(0);lc/t) ' ' , p (t) 
h (t) = -H— 



where we haveVjlsed the fact r that 



o:8 



' ~ Re * 



P M, V V 



Pr(t) 



1/3 



M (0) 



(2.3.1^) 



(2.3.19) 



\ 
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li (T) = 8.50ff7xlO"i - 1.7501xl6~ 3 3^+ 1.1142xl0~ 6 . ^ . (2.3.23) 

6 

and for the thermal conductivity of cthe fuel, we haye^ 

. kpCT) - 6.ll)2141 - 4.636522xl(f 3 1^+ 1. 306299xl(f 6 J*. ' . '(2.3,24) , 

The temperature is for the fuel in the case of bpt the moderator temperature 

otherwise. * . ■ . , 

Problem 2.3.^ . _* .' ' 

Given the following data: V — - • 

* ' ft 

- Rp = 0.015042 ft 

Fuel length = 12.00 ft . , - 

# * 4 Volume of reactor, vessel = 3643 ft - * B 

* 

Number of fuel assemblies = 145 
, Fuel rods per assembly 208 

• Rod pitch = 0.04733 



\ 



t p = 43.214 lbm/ft 3 



> 4 



'v = 16.3 ft/sec 
U -5.828 x 10" 5 *• 
D = Q. 04377 ft 
Re = 528,998 

» 

Pr = 1.01 

= 0.301Q Btu/hr-ftr ? F 



a) show the Cdlburn correction is a -0. 03148 

2 

b) h - 8232 Btu/ft -hr-°F. 
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v L 



X 



2.4 Numerical Solution of the Dynamics Equations 



I 



* < 



The numerical integration, of the differential equations can be carried out 
' once the various parameters, such as the heat transfer coefficient, thermal 

J 4 

conductivity, etc.., have been obtained. We again write the equations so that 
we can proceed to use Hansen's method' in their solution. FUMOTEM solves the 
• kinetics equations using Hansen's method.' The* system of equations is 



d£(t) 

— — = A i(t) + B 



\ 



(2.4.1) 



where 



i = 



p(t) 



Q(t) 



T p (t) 



^(t) 



B - 



0 
0 



V M l (t) 



p M V M 



(2.4.2) 



4 



and 



A = 



p(t)-e 

' A 



r 



2p F C F V F 



-X 



p F C F V F 



0 . 



- 4*1^ AtHc^L^ 



" P F C F V F 



p M V M 



(2.4.3) 
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subject to the constraint that 

P(t) = P Q 4: AT M + a F AT p . 



(2.4.4) 



9 



All the parameters involved in Equations (2.4.1) - (2.4.4) are defined in 
section 2.2. # % 

Hansen's method must be generalized slightly for the'solutiori of Equation 
(2.4.1). We break the A matrix into three parts: 



with 



D(t) 



P(t)-g 
A 



0 



L = 



3 
A 



>fVf 



P M C nV 



A = U+ D + y, 



0 . 



-A 



0 



P F C F V F 



ft) 



• 0 
0 

P M V M 

' 0 



(2.4-5) 
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and 



0 . > 0 



u = 



Equation (2.4.1) can then be written as] 



" § (t) ^ (t) = (L + U) i (t) + B. 



(2.4.6V 



We wish to develop an iteration procedure^ for the solution of this system%, 
so we begin the calculation at time t^ and advance to time t^, and define 



(2:4.7) 



Before we use "this relation though, let's multiply each term of Equation (2.4.6) 



by the integrating factor e o ■ 



-r D(t')dt' 



s'o , 



D(t:)dt' d± D(t*)dt' 

O = *T 0 = 

dt ■ e 



g(t)$(t) 



-r D0t') f dt' 



a+U)((.(t) + e J 



or 



A -f* D(t')dt* 
a , o = 



dt 



{e 



-/?D(t»)dt' ^ -/Jg(t')dt' 
.£(t)} = e (L + U)*(t) + e 



B. (2.4.8) 



The left side is an exact differential so we integrate from.O to tivand obtain 
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*-r fi(t')dt' . ih rh -r D(t')dt 



i'(t) 



h rn 
0 *°J 0 



,a+y)i(t)dt<+ / e 

0 



/ 



h -r D(t')dt' 



B(t)dt 



or 



' . +/".D(t , )dt' 
i(t 0 + h) - e • 



K /J D(t')dt' -/ C D(t')dt' • 

e e ° = (L + U)^(t o +9)d0 



• h D(t')dt* - /? D(t l )dt« 

0 = rt = 



+j e° 0 0 ,. B(t + 0)d0 . (2.4.9) 



Here 



and 



t o < 0 < t]L = t 0 + h 



d0 = dt. 



Notice in Equation (2.4.9) that 



/* D(t , )dt' - /' D(t')dt' '/J D(t')dt.' 



= e 



't+e S<t')dt' 

0 



Therefore the equation' that must be solved is 



r D(t')dt'~ 

■0 » 



i(t p + h) = e 



+ (t 0 ) + f e 



' h./J^DUOdf - 
' / e ° ... "(L + U) 



t +0 
o , 



i(t p + 0)dt 
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' t • 



.h /'; +0 g(t')df 

+ | e ° ' B(t Q +*0)dP . 

o 



(2.4.10) 



Equation (2?4.10)*is an integral ^equation so we must approximate its 
solution. To do this, we assume that , . 



0) 0 

i(t o + 0) - e 0 ±(t 0 * 



(2.4.11) 



where' 0)^ is'the largest eigenvalue of th^ matrix A evaluated at time t Q . This 
«. j. * » , * 

•means we must so^ye the equation 



A - o)I 



= 0 



or 



det 



P(t 0 )-B 



- - 0) 



2 P F C F V 



-A-0) 



p F c f v f 



o- 



>f c fV 



P M V M * 



= o, 



to obtain* the eigenvalues, to- , a) 1 u> 0 and a)_, where 




£2.4.12) 



P 1 2 3 
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The o) f s are, of course, time dependent quantities. ^ 

' We also look at J&(t + 0) and assume that we cah^expand it in a Taylor series 

* * 0 * *. Ki 

so.flhat 



' dB (t Q )' 2 d 2 B(t ) 

B(t Q + 0) = B(t o ) + 0 Tt ■ + f0 -^- + 



(2.4.13) 



and also 



i(t + 0) + U(t +0) = L(t ) + U(t ) + ... , (2.4.14) 

= 0 ■ 0 =0 = 0 



We shall keep only the 1st terms of jthese -^expansions . ^ 

Inserting Equations (2.4.14), (2.4.1J) and (2.4.11) into Equation (2.4.10) 

we get. (after assuming that /^D(t')dt' = D(£ n ) h ) ' 



6h' ; '. 



D(h-0) 



i(t Q + h) = e~-i(t 0 ) +J e? ■ (L + U) e u £(t ft ) d<3 



0) 0 

o 



•h g(h-e) ' 

+ / e Bft )d0 



— o 



Dh Dh 
•- e" £(t ) + e~-(u> I - D) 



(o> o ! - D)0 



,(L + U)j.(t Q ) ' 



' Dh — -D0 

K + e" (-g-) e = 



Putting the limits into this equation and simplifying we finally have 



N 



Dh , u hi Dh 

Ht + h) e" i(t ) + (w I - D)" i [e * %e~ llL^) + 2(t Q ) lAC^) * 



-1 Dh- 
+ D (e = - I)B(t 0 ) 



(2.4.15) 



If we now set 



and 



Then Equation (2.4.15) becomes 



(2.4.16) 



where 



Dh - -1 a) Ih Dh 

H = e + (w I -. D) [e ° = - e = ] (L + y) 



(2.4.17) 



with all quantities Evaluated at time t and 



-1 Dh 

R = D (e~" - I) .' . • 



(2.4. 4 18) 



We" new write & and B. explicitly for our problem. To do 'this, consider 
Equation (2.4.17) and the explicit relation for L and U; then 



H 



H l 


H 2 


0 


0 






H 4 ' 

/ 


Q 


Q 






l 

0 


H 6 


V 7 


9 




0 . 


0 


H i 
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where 



• H l 



P(t 0 )-B 



r '■ 



'h 2 =, k 



_• - P (t 0 )-e „• 

0) h : h 

e * - e * 



P(t 0 )-e 



H, 



6 -r w o h _ e -xtri 



-Ah 



I. "5 



4tt 



o) h 
o 



- e 



p F C F V F 



. o 



4tt 



p f c f v f 



H, 



47rk' F L^ 

#7 



H. 



p F C F V F ' 



oi h 
o 



P F?F V F 



0 ~ p F C F V F 



0) - 
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a) h 
o 



2ft 



-t e 



"mYm 



+ 2 *M 



2n 



p M V M 



MM J 



u) h * D V' 
a o" ■ M ,M 
e - e % 



0) + 



2i Vi 



0 P M V M 



and the R matrix is 




1/1 " Ah N 



" . _4irkpLpNh - 



^M 



-2 



(1-e 
{ 



% V M 



(2.4.20) 



The basic iteration procedure in J5UM0TEM is as follows: 
1. Calculate the pertinent reacto^ parameters such as a^, h^» mr, P^, P, C^, 
k^, and for the guessed initial conditions. The parameters read in are 
^Ml* Ct M* k> ^ 0 > v (the coolant* velocity) , some" optional ^settings for print 
out and 'the t-ime length tfie program is to operate. 
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Construct the vector $_(0) . 'Ve can choose" the parameters of interest, i.e% 



P(0) = 500 MW, . . • y 

and i 



Q(0) = P(0) MW. 



The fuel temperature T F (0) and the moderator temperature ^(0) are completely 
determined by P,(0) and the above read-in parameters. Tp(O) and ^(0) are 
difficult to obtain since C^, k^, y etc. depend on the temperature. This 
difficulty is overcome by using an iteration procedure. The method is ah # 
follows: v 

a} From P(0) and initially guessed values t£ 0) (0) and T.^ (OX 



calculate K, from the relation 



\ 



v 4 ^ Equation' (2.3.20) is used to obtain P^C^) from th - e g uessed 

(0} * I ' 

moderator temperature T _/ ,(0) * Equation (2.2.15) is now 

used to obtain *an improved guessed power ?^ (0) 

P (1) (0) = 2^' (T ( ° } (0) - T M1 ). (2.4.22) 



b) The fractional difference between *tfie actual power and P^j[0) is 

p = I p (1) (Q» - p(o) L ./ 



P(Q) 



■^4.23) 



If E <_ A (A= 0.01 aad is an input) then the guessed temp'erature 
T^(0) is the correct moderator temperature. ,If * 
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E '> A (2. A. 24) 



then change (0) to T^(0) and repeat the above until 
convergence is achieved . ; 
•c) The T^(0) is now, use <%in Equatffn "(2.2.17) which*, for equilib- 

> * « 

rium, becomes ,* after solving for T_,(0), 

n 

V°'):V 0) + -8^ P(0) - * 02.4.25) 

This step also requires an iterative procedure since 
is a function of the f\iel temperature as seen from Equation' T2.4.25) 
Determine the largest eigenvalue of the equation 

|A - oil] f 0; 



This will be* thte solution of a 4 x 4 determinant which is rather easy on 

the computer. From this determine the largest root u) . The Newton-Raphs"on 

method is used, to calculate all four roots of this polynomial and then' the . 

largest root is picked by comparison of the rootg. The subroutine POLRT 

is usted for the determination of u) for each time step. Module RD-1 

o 

describes' the Newton-Raphson method. 

Construct the>H(t ) matrix using Equation (2.4J.19) and R(t ) usiing (2.4.20) 
m »o • n ~ o 

Determine the vector^ from Equation (2.4.16), i.e. . ** 



in -'H(t )jL + R (t > B 



Repeat the above procedure using step 3 and continue as long as required to 
achieve ^the^oliition over the time domain of interest. ^The -time steps are 
chosen by the criterion h -inland h is constrained to he in the Interval 

0 a f « 

0.0Q5 < h ^ 0.05 sec. ' » ' 
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2,5 Input-Output Data" for Code FUMOTEM 



The -Input data required for the program are presented below: 



Data Card Format 
Number Statement 
Number 



Format , Unit Variable 
Name * 



1 


20 . 


12' 




N0PLT 


2 

* 


30 

-> 

* 


• 

F10.5 
F10.5 
F10.7 
F10.7 
110 


~"-l 
gec_ 1 

sec 
$ - 


BETA 

•x 

XL 
RO 
NRO 






F10.5 


•sec 


RTIME 






F10.5 


-1 

sec 


A 


J 




pin A 


o F "l. 
r 


AM 

Ail 






FIO.6^ 
F10.6 




G 

PTFO 


4 - 


50 


F10.5 
. F10.5 


ft 
ft 


RF 
PC 






F10.5 

F10.2 
. F10.2 

no 

110 


Itu 

lb-°F, 
lb/ ft 

ft 


CPFX 

T7TM7XTC 

FH 
NA 
NFRPA 


5 • 


■ 60- 


F10.2 
F10. 2 


ft/sec 
°F 


v" 

• TMI 


6 


' 70 


Fl.0.4 - 
F10.4 
' F10.4 


sec 
sec 


TE 
DH 
DELTA 

♦ 


7 


... 80 


F10.2 
F10.2 


MW 
°F 


PPW 
TGUES 






' F10.Z 


°F 


TGF 



Description 

Plotting -option 

1 — plot " , 

0 — no plot 

Delayed neutron fraction 

Delayed neutron decay constant 

Neutron generation, time 

reactivity, p Q 

option for the type of ^ 

reactivity insertion 

Time duration reactivity is 

inserted for ramp input (NR0=2) 

reactivity insertion rate 

Moderator temperature ^ 
coefficient of reactivity » l 
Constant for ^resonance capture 
Resonance escape probability, P 



Fuel radivjs, Rp 

Distance between fuel pin 

centers, P 

Specific heat of the fuel, 



Fuel density, p p 
FueJ. height, L 
Number of fuel assemblies 
Number of fuel rods per 
assembly 

Average velocity of coolant, v 
Inlet coolant tempetature to the^ 
core T M1 . - 

End of calculation time 
Time increment 

Convergence criterion for the 
calculation of the initial - 
equilibrium state 

Initial equilibrium power 

Guessed initial moderator temperatute 

to obtain equilibrium 

Guessed initial fuel temperature for 

obtaining -equlibruja conditions* 



f 




' 0 0 0 0 1 1 1 6J 0 0 0 0 0 0 1 t 0l 0 0 0 0 0 0 l I 0l 0 0 0 0 0 0 0 ttOOO 00 OpO 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 O'OO 0 0 0 0 0 0 0 0 0 0 0 0 0 

I I } 4 S t 7 1 * 10 II II IJ M 15 II 17 II II 30 31 33 31 34 25 3S 27 21 31 70 II U 11 )4 15 JO JMI IS 40 41 43 41 44 45 41 47 41 41 50 SI 53 51 54 S5 5407 31 SI 10 II II IS M 13 II 17 $1 II 70 71 72 73 74 75 71 77 71 7J 00 . 

1 1 1M 1 1 jl 1 1 1 1 1 U 11 i 1 1 1 1 1 1 11 1 Vl 11 1 t»1 U If 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 H 1 1 1 1 11 1 1 ' ' 1 ' ' " 
5 2 2 22-2 2 222 222 2 23 2 2 2 2 2.22 2 2-2 2 2 2 2 2 2 2 2-2" 1 2 222 2 2 2 2 222 2*2222 2222 2 22 2 2 2 2222 2 22 2 222222 222 2 2.V 
I 3 3 3 3 3 3 34 3 3 3 3 3 3 3 3 3« 3 3 3 3 3 3 3 3 3 » 3 3 3 '3-3*9 3 3 3 3 3 JS 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 f 3 3T3 3 3 3 3 3 3 3 3 3 3 3 3 3 ^ 
i 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 -4 4 4 4 4 A 4 4 4 4 4 * 4 4 4 4 4 Wijftt* 4 4 4 4 '4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 '4 4 4 4 4 4 \ 

a. ( - *.Vtt • * ' - 

8,55555 55 55 5 555555555 555 5 5.S 5 5 5 55 55 55 5 55S55ljisjSS555S55555 5555 5 555 5 5 555 5 5 55 5555 5 5 5 5 * 
| 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6.6^ 6 6 6 6 6 6 6 6 %\ 6 %% 8 6 8 6 6 6 8 6 6 6 6 6 6 6 6 6 6 6 6 S G 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 . 
| 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7-7 7 7 7 7 7 7 7 7 7 7 3 7 1 lSl1 7 T ; 7 7 7 7 ? 7 7,7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 1 7 ,3 | 
| it 1 8 8 8 8 9 8 8 8 8 8 8 8 8 8 9 ? 8 8 8 6 8 8 8 8 9 8 8 8 8 8 8 8 8 1 8 8 8 8 8 8 8 J 8! 8 8 8 bHh 8 8 8 8 8 8 8 8 8 8 8 8 8 8 SlY>HSSS 8 1 | 
9 99 91999 9 9 9 999 9 9 9 999999 99 9 9 9? 9 9 9999 9999 999999 999999 9999 99999 9 9999999 9999 9999 999 9 



In addition to these seven data cards, there is an eighth card which has 
-the impression as shdwn below. This 'card is necessary for the plot routine. , 

/* fH FRP1LM V: • r *~ T " ^ : ' ^ 



pB2 9361 



4- 
4i 



s 

. o 
«, o 



!!!!!! 

1111 IH 1 1 1 1 1 1 7 1 1 1 1 1 1 1 1 1 1 1 mi 1111111111 111111111111111 lllfl 1111111 1111111111111 
'2 2 2'2 2 2 2 2 2 2 2 2 2 2/2 2 2 2 2*2 2 2 2 2 2 2 2^ 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 22 2 2 2 2 2 2 2 2 2 2 2 * 

33333333.3331*3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 _33 * 
4^4 4 4 4 4 4 4 . 4 4' 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 ^4^4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 44 4 44 j 

5 5*.5 5 5 5 5 5 5 5 5 I 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 55* 5 5 5 W*Ji 5 5,5 S S 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5555JV555 5 

6 6 66. 6 6615 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 56 66 6866>{6 8,86 8 1 6 6 6 6 6 6 6 6 6 6 6 S 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 8 S i 6 B - 8 
-.77 7 7 7.7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 77 7 7^-71 7 77*7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 777 1\V> \ 

_ J88888888888 8 8 888888688888888 8 8 88888888888888888 888888888 6 8 8 8 8 8 8 8. | 
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The output of 9 FUMOTEM consists of three parts. The first is-simply a 
writing of all the input data. £he second is the equilibrium state calculation 

• r • * , 

• and the third block of data is the fuel temperature, moderator temperature • 

reactivity, exit temperature, power, a p and Q as a function of time. If the 

* 

J user specifies, a plot of power , precursor power, 1^, T M and reactivity as a^ 
function ftf time i§ provided. . , 



Problem 2". 5,1 

'ftun FUMOTEM for the sample data cards shown. 

Problem 2.5.2 

Run FUMOTEM for the reactivity dnput 



p(t) - $0.50 cos 2,5t 



with the following parameters: ' 

= - 0.00005 



.NOPLT m 1 

V = 0.00645 

\ m 0.07695 sec 

\ 

RTIME - A.O sec 



-1 



y = 0.0002 



P =0.80 
o 



P 



0.0150A ft 
0.0A733 ft 



C fc = 0.0590 Btu/lb-°F 
P F •« A3. 2 lb/4t 3 
L ■• 12 ft 

t 

TE « A.O sec 



Number of fuel assemblies - 7A5 



Number of fuel pins per assembly - 208 



T M1 " A0 °° R 

i 

DELTA « 0.01 



DH = 0.01 sec 

P = 1000 MW 
o 

TGUES.= 200. 0°F 
TGF - 500. 0°F • 

= 13.0 ft/sec 



( 



Also run the same calculation for v = 26.0 ft/sec. a 

c 

FUMOTEM is written in FORTRAN in single precision except for the eigenvalue 
calculation. The solution of the equation 

ia - a) i r» o 

for the root w Q , u^, u> 2> is done in double precision. All four roots are 
determined and* the largest one is picked to form the H and matrices.^ 

The memory requirement is about 40 kilobytes and the execution time varies 
with the time length. X^he reactor is to be simulated. Generally it takes 
about 1 1/2 seconds of computer time to simulate one second of reactor 
transient time. 

• ; • / 



U2 
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program. 



* * * . List of Symbols for FUMOTEM 

The following symbols are listed in alphabetical order in the FUMOTEM 



/ 



A 

AF 

AM 

AMTRX 

AMX 

AX 



BE 

BETA 
CK 
t 

COL 
CP 

CPF 
CPFX 

CPG, CPMX 
CX, CGUES 
DE' 

DELTA 



C 

C PM 
HC PF 

<W 

*M 
De 

A 



Period of reactivity insertion 

Fuel temperature coefficient 

Moderator temperature coefficient 

Subroutine to form the A matrix 

Elements of the A matrix 
iT 

Coefficients of *the EIGEN4 polynomial, 

.5 

c ' n 
i.e. , if ) a a) , the a . 

* • 

B matrix ^Lement, 

Delayed^ neutron fraction 

Function subroutine to calculate the 
water thermal conductivity as a 
function of temperature (Equation 
2.3.21) 

Colburn number (Equa tijpn 2.3.17) 

i 

Heat capacity of moderator as a 
function of T^ % * 

Total heat capacity of the fuel 

Heat capacity' of a single fuel rod * 

Total heat capacity of the water 

Thermal conductivity of theyater ^ 

Equivalent diameter (Equation 3.3 T .16) ' 

Convergence criterion for equilibrium 
calculation 
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DF 

DH 

DH1 

DH2 

DM 

DOTMG . 
DX, DGUES 
EIGEN4 

E-IGENV 
FA 

FCA . 

FDENS 

FH 



FK 



FLA . 

FX, FGliES 
G 

GXN 

H 
HH 

HMTRX 



2R F 



At 



t 

P M 

to 



•••)»* <^ 
15k 



HGUES 



Y 



H 



• h T 



Fuel diameter 

Time increment (TINCR) 

Lower limit time increment 

Upper limit time increment 

Time at any instant 

Mass flow of the coolant 

Density of the coolant 

Subroutine to change, | A - *o)l| == 0 
to polynomial form 



\ 



S 



Eigenvalues of the A Hfatrix 

. i 

Heat transfer area 
Fuel cross section area' 
"Fuel density ^\ 
Fuel rod length 

Function subroutine used to calculate 
\ the fuel conductivity as a function ^ 
of temperature 

/ | Flow area 

Thermal conductivity bf fuel 

Constant in Equation (2*3.12) 

Subroutine) to multiply an nxn matrix 
with a column matrix 

Heat transfer coefficient 

Elements of the H matrix (Equation 2.4 

Subrou^ifie to form the.H matrix 

Element of the column matrix H<|> 

fi ; 

Heat;%r£n8fet coefficient 



MTRX 




MX 


t 


N 




NA 




NCHAti 




NFRPA 




NIN 




NN 




NOPLT 


• - 






NRO 


« 


NW i 




NZ * 




• 

PAREA 












pi 


TP" 


PKL 




POW 


p 




\J 


"D-D * 
PP 


r 


' PR 


Pr 


*PW 


P 




w 


QQ, QW 




■ RIN J . 





. Order of the'nxn matrix 



RB V 



Dimensioned time 

Number of iterations 

Number of- fuel assemblies 

; Total number -of coolant channefs ' 

Number of fuel rods per assembly 

Number of iterations calculated 
frop insertion time- and initial' 
time increment 

Total number of iterations, calculated 
from end time and initial time ^increment 

Option for plotting ' ' ~\ 

READ statement unit number 

Option for reactivity insertion » 

WRITE statement unit number 

Dummy variable used for print out 

> 

Square pitch area £or channel^ 
distance between fuel pins 
3.14159 t 

4lrk P L F 

Input^pbwer in MW 
Sower in Btu/hr • 

Prandtl number function ^subroutine 
Wetted perimeter < 
Precursor density (Power) ^ 
Reactivity inserted • 
Element of column matrix R B 
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RCOS 
RE ' , ( 

RHO 
RL 

RMTRX - 

RN • 
RO 

RSIN 

RP 
RR 

RTIME 

RX - 
RY 

TCFA 

TCFLA 

TE r 

TFA 

TFG 

TPO 

TGUES, .'TMG 



p cos a t 
o 



Re 



p = P Q (l+a t) 



p = p sin a t 



Pr 



P in + P i 



r . 



T 



Function subroutine to calculate 
reactivity for a cosine insertion 

Function subroutine to calculate 
the Reynolds number 

Fuel radius k , 

Function subroutine to calculate 
density of H 2 0 

Function subroutine to calculate 
reactivity for a ramp 

Subroutine to form the R matrix 
(Equation 2.4.20) ^ 

Reynolds number * 

Reactivity 4-nsertad at 4: ^0 

Function subroutine to calculate 
reactivity for a sine insertion 

Prandtl ^umber 

Total reactivity, inserted plus 
feedback 

Time when inserted reactivity is 
removed (for rltrop reactivity only) 

, Element of R matrix (Equation 2.4. 20) 

9 

„ 4 

Dimensioned reactivjLty 

Tota^ cross sectional area of ^ the fuel 

* Total cross Stecfcional area of coolant 

End of calculation time 

Total heat transfer area 

Initial guess for fuel temperature 

Equilibrium fuel temperature corresponding 
to P % 



•Initial 'guess for moderator temperature 



0 



TMI * T * * Inlet coolant temperature 

TMO - T * * Equilibrium coolant temperature 



corresponding' to a given P 



o 



TMOUT ' Exit temperature of core water 

TPLOT z Subroutine to plot five variables at . 

the same time with respect to the 

* " independent variable time 

» - " 

U $ e \i Function subroutine to calculate t^e 

^ - * viscosity of water as a function of 

* temperature 

UGUES, UX .» Viscosity, of 

V ( V f Mean' velocity of coolant 

* * 

VF V * ^ ' - Fuel volume 

• . F 

w 

VM e > • V Coolant volume 

* ( 

w The latgest eigenvalue to the 

; o equation. |A - a>l| = J ' 

X X . Decay constants of the' delayed neutron 

* - ' 4 group 

vt A The neutron generation time- 
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Flow Chart for FUMOTEM 



^ - Start ^ 



s 

V 

( • 



0 




Geometry 
Calculation 



i 



Equilibrium 
conditions 



Form Initial <J> 



Store 6 for 
plot 





Calculation to obtain 
^'fuel area, flow' area 
etc. are necessary 



Once P is read in, 
then* a consistent fuel 
and moderator equilibrium 
temperature is obtained 



Initial condition 
vector is fon^fl from 
the equilibrium values 



r 



0 



r\ 



•5- 



ii9 ■, 



* 



42 
t 




If the feedback reactivity 
coefficients are positive, 
the program stops* 
Actually, cy* 9^ must be 
less than or equal to zero 



Po (t)=P o 



T 



P 0 (t)=p 0 (l+at) 



-0 



■ © 



0 (t)=p' cos at 
K o o 1 



± 



7777 



P 0 (t)=p Q sin at. 



One of these 
four reactivity 
insertions is 
to be used'' - 



Form B 0 



Form A 

rs 

matrix 



r 



Really, since the B 
<vector has 3 zeroes, 
only the number 

2 \ t hi (t) 

needs to be solved. 



EIGEN4 # 



This calculates the 
I A - wl| = 0' and puta_^ 
it into polynomial 
form 




Calculates the largest 
eigenvalue w 
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© 



t 


Add 
Matrix 
to obtain <|> 


J 








J 


P(t) - P Q 


(t)+ p f (t) , 







At is constrained 
to remain in the 
▼ time interval 

0.005 < At < 0.05 1 



GXN is a subroutine 
which multiplies 
mattices to obtain 
* (t) 




Calculate total 
reactivity since 
knVw inserted and 
* feedback components/ 



If t > t , set 
8 V'P(t)- = Pf(0 



If not at ^jid of 
calculation, go 
back ' 



12\ 



20 



30 



//WATflV 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c • 
c, 
c 
c 
c 
c 
c 
c. 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

c 
c 

-c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

t 
c 
c 
c 
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0NEGA,TlM£»300tPAGES*50 ' 
PWR FEEDBACK 



C00E NAME 
OBJECTIVE 



FUMOTEH . 



PROGRAM 



REACTIVITY . 

WRITTEN IN SINGLE PRECISION . 

DESCRIPTIONS OF INPUT PARAMETERS . 
FUNCTION 



FORMAT- PARAMETER 
NUMBER 



NOPLT 



BETA 



X 

XL 
RO 

NRO 



\ 



40 



RTIME 
A 

AM 

G 

PTFO 



OPTION FOR PLOTTING 
1 - PLOT THE RESULT . 

0 - NO PLpt . 

AVERAGE FRACTION OF DELAYED NEUTRONS 
IS 0.00645 t ONLY ONE DELAYED 
NEUTRON GROUP IS CONSIDERED • 
DELAYED NEUTRON DECAY CONSTANT . 
NEUTRON GENERATION TIME . 
INITIAL REACTIVITY INSERTED . 
INPUT RO IN DOLLAR UNIT . 
OPTION FOR REACTIVITY INSERTION 
AS A FUNCTION OF TIME . 

1 - CONSTANT REACTIVITY . 

2 - LINEAR RAMP- INSERTION • 

3 - COSINE VARIATION OF REACTIVITY . 

4 - SINE VARIATION OF REACTIVITY 

WITH TIME 
INSERTION TIME 
PERIOD 9 CONSTANT . * 

• > 
MODERATOR TEMPERATURE COEFFICIENT 
CONSTANT GAMMA IN EQUATION 2.3.13 
OF RD-2 

RESONANCE ESC/fPE PROBABILITY . 



FUEL ROD PROPERTIES • 



50 



j 



60* 



RF 4 RADIUS • 

PC DISTANCE BETWEEN RODS • 

'CPFX SPECIFIC HEAT . 

FOfcNS DENSITY . 

FH LENGTH ^ 

NA NUMBER OF FUEL ASSEMBLIES • 

NFRPA fUEl* ROPS PER ASSEMBLY . 

V ' VELOCITY OF MODERATOR/COOLANT 

TMl INLET TEMPERATURE OF MODERATOR 



4 


A 


i 
i 




A 
A 


o 
C 




A 
A 


* 3 




A 

A 


A 




A 


c 
J 




A 
A 


A 
O 


a% ft J r% A lift 

PWR AND 


A 


( 




A 

A 


Q 
O 


TURE 


A 

A 


o 


ON OF 


A 
A 






A 

A 


1 I 




A 

A 


1 O 




A 

A 


13 




A 


1 A 




A 

A 


1 c 
1 5 


• 


A 

A 


1 A 




A 


17 




A 


18 


UNIT 


A 


19 




A 


o n 
Zu 




A 


O 1 

Z I 




A 


o o 
ZZ 




A 


23 




A 


24 




A 


o c 
Z D 


i 


A 

A 


OA 

do 




A 
A 






A 
A 


OQ 

Zo 


1/SEC • 


A 
A 


OO 


SEC 


A 
A 






A 

A 


3 1 




A 

A 


* o 
oZ 




A 


33 




A 


34 




A 

A 


3D 




A 

A 


3d 


* 


A 

A 


3 f 




A 
A 


30 




A 
A 


3 7 


b cC 


A 
A 


*tv 


1 /b CL 


A 
A 


* I 


i 

• 


A 
A 


AO 
**Z 


1/F« 


A 
A 


L.1 

TJ 




A 
A 


AA 


c f\&T ret 


A 
A 






A 
A 


AA 
•to 




A 


47 




A 


,48 


•* 


A 


49 


FT • • 


i A 


50 


FT. . 


A 


51 


BTU/LB-F 


A 


52 


LB/.FT**3 


A 


53 


' FT . 


A 


54 




A 


55 


• 


A 


56 




A 




ft/sec. 


A 


58 


F. 


A 


1 59 



I 



123 



46 



7 
8 
9 
10 
li 



12 
*3 



C 
C 
C 
t 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

' c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

, c 
c 
c 
c 
c 
c 
c 
c 

c 



c 
c 
c 
c 



70* 



80 



TE 

OH 

DELTA 

PPW 
TGUES 

TFG 



END TIME OF CALCULATION* 
TIME INCREMENT t THE RANGE MILL 
BE BETWEEN 0.005 TO 0.05 SECOND 
CONVERGING FACTOR . 



SEC' 



GI VEN 'EQUI LIBRIUM POWER AT T=0.0 MW 
COOL ANT/ MODERATOR GUESSED TEMPERATURE. F. 
CORRESPONDING TO PP . 

FUEL GUESSED TEMPERATURE pORRESPDN - F. 
DING TO PP • ' ^ 



MAIN PARAMETERS OTHER THAN INPUT . ^ a f 



PDWER. DENSITY- . 
PRECURSOR DENSITY • 
FUEL TEMPERATURE • 
MODERATOR/COOLANT TEMPERATURE . 
REACTIVITY AT ANY INSTANT.. 
TIME ELAPSED • 

TOTAL NUMBER, OF ITERATIONS • 
COOLANT/MODERATOR EQUILIBRIUM, 

TEMPERATURE . 
FUEL EQUILIBRIUM TEMPERATURE . 

HEAT TRANSFER tOEFFICJEJiT . 
TDfAL MASS FLOW 

DUMMY INDICATOR FOR PRINT-OUT THE 
EQUILIBRIUM CONDITION BEFORE AND 
AFTER REACTIVITY INSERTION °. 

FUEL 'TEMPERATURE COEFFICIENT . 

NUMBER OF ITERATION AT ANY TIME 



PHI (Dt MPEt MPEL 

PHI (21 tMOU 

PH1(3),MTEEF 

PHI (4)tMTEEM' 

RR,*Y 

DMHf MX 

NN 

TMO 
TFO 

HXf HGUES 
DOTMG 



HI 

AAF 

N 



REAL MPEL(5OO),MPE(56o),MQU(5OO)tMTEEF(5OCI,MTEEM(5OO),MX(5O0)',RY( 

1500) • , 
DIMENSION* AMX (4 t^lfHHl4«4> f RX<4» 4) 
DI MENSI OH R&OTI (4) , EI GNV( 4 ) ,COF < 5) t AX ( 5 > 
DIMENSION PHI (4) , BE ( 4 >.,R6 ( 4 1 f HP< 4 ) 

COMMON BET AfXLtXt FX? RFfOHt FDENStCPFt VFtFHt VM f DXtppTMGt WOtHXtRR tC 
lPMX»PItNROD , 
DOUBLE PRECISION AX »COF, E IGNV ,ROOTI 



l define reynolds number t prandtl number , heat transfer co 

Efficients , fuel coefficient temperature and reactivity . 

RE (RElf RE2?RE3?RE4)=REl*RE2*RE3/RE4 
PR(PRl*PR2t PR3l*PRl*PR2/PR3 

H(Hl-,H2tH3tH4,H5)=Hl*H2*(H3«*l8)«(H4**.3333)/H5 
» AF(AFl»AF2»AF3)=-AFl*AL0G<1.0/AF3)/SQRT(AF2) 
RT I RTl t RT2 tRT3f KT4t RT5? RT6t RT 7 )*RT 1+ ( RT2*RT3)*RT4-RT5*RT2*R T6+RT7 



—FUEL ANb WATER PROPERTY AS TEMPERATURE DEPENDENT . 



FK(T)«0.61021E1-0.46365E-2»T+0.13063E~5*T**2 
CK<n*0.1l7lU0«l39l0E-2*T-0.18102E-5*T**2 



A 


60 


A 


61 


A 


62 


A 


63 


A 


64 


A 


65 ' 


A 


66 


A 


67 


A 


68 


A 


69 


A 


70 


A. 


71 


A 


72 


A 


73 


A 


74 


A 


75 


A, 


76 


A 


77 


A % 


78 


A* 


79 


A 


' 80 


A 


81 


A 


82t 


A 


83 


A 


84 


A 


85 


A 


86 


A 


87 


A 


88 


A 


89 


A 


90 


A 


91 


✓ A 


92 


A 


,93 


A 


94 


A 


95 


A 


96 


A 


97 


A, 


98 


A 


99 


A 


100 


A 


101 


A 


1D2 


A* 


103 


A 


104 


A 


105 


A 


106 


A 


107 


A 


108 


A 


109 


A 


110 


A 


111 


A 


112 


A 


ILi. 


* A 


114 


A 


115 


A 


116 


A 


117 


A 


118 


A 


119 



124 



47 



14 
15 
16 



17 
18 
19 

20 
21 
22 
23 
24 
25 
26 
27 
28 
29 



30 
31* 
32 
33 
.34 



RH0(T)*0.57788E2*0.28018E-1*T-0.883.46E-4*T**2 
UtT) =0.85067-0. 17501E-2*T*0. 11420E-5*T**2 
CPIU«0.47088E1-0.15753E-1*T*0.17233E-4*T**2 



35 



36 
37 



38 
39 
40 
41 
42 
43 
44 
45 
*r6 
47 
48 
49 
50 



C 

C 

c 
c 
c 



c 
c 

. c 



c 
c 
c 
c 
c 
c 

c 
c 

-c 
c 



c 
c 



- MAKE SURE THAT NN EQUAL TO DIMENSION NUMBER OF 
THE PLOTTED VARIABLE . 



NN=1600 
NR=5 , • 
NW*6 * 
NZ=0 

PI=3. 14159 
N=l 

DM=0.0 
*DH1=0.005 
DH2=0.05 
0MH=0M*3600.0 
RaO.O 
RR-0.0 
MX(1)=0.0 



-#-- READ INPUT DATA 



READ (NR 1 29) NOPLT * 
READ (NR»30) BETA, X ,XL , RO,NRO ,RTI MF, A 
READ <NR,31) AM,G»PTFO 

READ (NR,32) RF , PC ,CPFX t FOENS ,FH*NA, NFRPA 
READ,CNR|33) V,TMI 



THE INITIAL VALUE OF OH WILL CHANGE ACCOROING TO THE LAR 

GEST EIGEN VALUE OF A-MA'f R IX . 



READ (NR t 34) TE »DHtDELT A 



READ INITIAL POWER Of SIRED* AND GUESSED FUEL AND COOLANT 

TEMPERATURE . 

READ <NR,35) PPW, TGUES, TFG 
PP»3.412D6*PPW < \ 



-PRINT OUT INPUT DATA AND THfc INITIAL . 



WRITE (NW,36) 
WRITE (NW»37) 

WRITE <NW,38> BETA.X,XL,RO 
R0«R0*8ET A 

WRITE (NW,39) RTIMEtA 

WRITE (NW,40) AMtGt PTFO 

WRITE INW,41> RF,PC,CPFX,>=OENS»FH 

WRITE*(NW»42) NA»*NFRPA % 

WRITE (NW f 43) V,TMI 

WRITE (NW«44) NOPLTtNRO 

WRITE INW.45) PPW,TFG , TGUES 

WRITE (NW.46) TEfDH 

WRITE (NW»47) DELTA 

CHANGE SECOND TO UNIT HOUR 
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121 


A 
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A 


123 


A 


124 


A 
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A 
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A 
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A 
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130 


A 
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A 


L32 


A 


133 


A 


134 


A 


135 


A 


136 


A 


137 


A 


138 


£ 


139 


A- 


140 


A 


141 


A 


142 


A 


143 


A 


144 


A 


145 


A 


146 


A 


147 


A* 


L48 


A 


149 


A 


150 


A 


151 


A 


152 


A 

« 


1 51 


A 


1 54 


A 


155 


A 


156 


A 


1 57 


A 


158 


A 


159 


A 


160 


A 


161 


A 


1 62 


A 


163 


A 


164 


A 


165 


A 


l£6 


A 


167 


A 


168 


A 


169 


A 


170 


A 


171 


A 


172 


* 


173 


A' 


174 


A 


175 


A 


176 


A 


177 


A 


178 


A 


179 



125 



- 48 - 



C - 



51 
52 
53 
54 
55 
5$ 
57 
58 
59 
60 
61 



62 
63 
64* 
65 
66 
67 
68 
69 
70. 
71 
72 
73 
74 
75 
76 
77 



78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 

^90 
91 

92 
'93 
94 



95 
96 
97 



C 

c 
c 
c 

*2 



C 
C 
C 
C 
C 

3 
4 



N0H=*IFlXJTE/0H*DH/2.0) 

IF (NDH.J_E.NN) GO TL 1 

WRITE <N)f, 48) 

GO TO 28* 

XL*XL/3600.0 

X=X*3600. 

V*V*3600. 

DH1=DH1/3600.0 

OH2=DH2/3600.0 • 

DH=DH/3600. 

A=A*360Q>0 



-GEOMETRY CALCULATIONS FOR SQUARE PITCH 



OF.*2.0*RF 

FA*PI*DF*FH 

FCA=PI*DF**2/4.0 

PAREA=PC*PC 

FLA*PAREA-FCA e . 

NROO*NA*NFRPA 

NCHAN^NROD * - 

TFA=NROO*FA 

TCFA=NRQD*FCA 

TCFL A=NCHAN*FLA 

VF=FCA*F**NROD 

VM=TCFtA*FH 

PW=PI*DF? v o 

DE*4.0*FLA/PH 

CPF'CPFX*FOENS*VF 

COL=0.042*PC/OF-0.024 

* ITERATE FUEL AND/OOOL ANT TEMPERATURfc UNTIL IT CONVERGES Tu 

THE CORRESPONDING P©WE|t AND ITS PRECURSOR . 

CONTINUE- \ 

UGUES=U<TGUES) j 

OGUES*RHO(TGUESJ \ 
OOTMG=DGUtS*TCFLA*V 1 
CGUES=CK(TGUES) 
FGUES=f;K( TFG) 
PKL=4.0*PI*FGUES*FH 
CPG=CP< TGUES) / 
TMG s TMl «-PP/ < 2»0*CPG*DOTM6) 
TFG=TMG*PP/ <2.0*PKL*NR0D) 
RN»R£(0G0EStVt0Ett*6UESI 
RP=PR(CPG»0(*UES,CGUES) ) 

HGUES=H(COL,CGUES,RN,RP,DE) • - A 
PGUES*2. D*CPG*OOTMG*ABS ( TGUES-TMI) 
OP=PGUES-PP 
OTM*TGUE$-TMG 
OELPG*ABS(OP>/PP 



SET THb GUESSED POWER ANO ITERATE UNTIL CONVERGE' TO THE 

CORRESPONDING FUEL AND MODERATOR TEMPERATURE ACCURATE , fQ THE 
VALUf-OF DELTA ^ ' " 

IF (DELPG-DELTA) 6,6,3 

IF (DTM) .4,6,5 . , ' 

TGUE S*TGUES*ABS ( DTM) /2. 0 



180 

181 

182 

183 

184 

185 

186 

187 , 

188 

189 

190 

191 

192 

193 

194 

195 

1.96- 

197 

198 

199 

200 

201 

202 

203 

204 

205 

206 

207 

208 

209' 

210 

111 

212 

2li 

214 

215 

216 

217 

218 * 

219 

220 

221 

Z22 

223 

224 

225 

226 

227 

228 

229 

230 

231 

232 

23* 

234 

235 

236 

237 

238 

239 



126 



- 49 



98 
99 
100 
101 
102 
103 



104 
105 
106 
107 
108 



109 
110 
111 
112 
113 
114 
115 
116 
117 
118 



119 
120 
121 
122 
123 
124 
125 
126 
1-27 
126 
129 
130 
131 
132 
133 
134 
135 
136 
137 
138 
139 
140 



141 
142 



C 
C 

c 



c 

c , 
c 



c 
c 
c 
c 



9 

10 



C 

c 
c 
c 

11 



GO TO 2 

TGUE S-TGUES-ABS<0TM)/2.0 
GO TO 2 

CONTINUE t 

QQ*8ETA*PP/<X*Xl) 

QW=QQ/3.412E6 



— INITIALIZE MATRIJTPHI 



PHM1)*PP , 

PHI ( 2)=0Q 

PHH3)=»TFG 

PHH4)-TMG 

DOTMG*DGUES*TCFLA*V 



STATE 



SETUP THE FUEL AND COOLANT TEMPERATURE AT EQUILIBRIUM 



FUEL AND MO 



TMO 3 TMG 
TFO=TFG 
RS =RR/8E TA 
WRITE <6,49) 

WRITE <NW f 50) PPWtOHtTFOt THOiRS 
WRITE (NW»51) RN,RP,COLt DOXMGfHGUES 
WRITE (.NW,52) 
WRITE <NW,53) 

CONTINUE / 
MTRX5^— ' ^ 

STORE THE VALUE OF POWER t PRECURSOR DENSITY 

DERATOR TEMPERATURE • AND REACTIVITY *0R PLOTTING 

PPW=PHH D/3.412E6 , / 

QW=PHI(2)/3.412E6 
TFG=PHH3> 
TMG*»HI<4) 

IF INOPLT.CQ.O) GU TO 8 
MPE<N)»PPW 
MPEL(N)=AL0G10< PPW) 
MQU(N)=QW 
MTEEF < N) 3 ,TFG 
MT£EM<N)=TMG. 

A-Af^AF < G» PHU 3) t PtFO ) % 

IF' i AAF.,LT.O.O.AND.AM.LT.O.O)^ GO TO 10 

IF ( A AF+AM ) 10 f 9»28 

WRITE <NW,54) 

UX=U(PHU ? 4)) 

DX=RH0tPHI<4)> 

CX«CK(PHI(4)> 

CPMX*CP<PHU4n 

FX=FM PHU 3) ) 9 

DOTMG=DX*TCFLA*V 

HX^HICOL tCXt RE( DXt Vt DE t UX) » PR(CPMXtUXtCX) #DE) 
TMOUT=2*0*PHH4)-TMI 

COMPUTE THE' REACTIVITY INSERTED AS A FUNCTION OF TIME 

. % AND REACTIVITY . ' 

GO TO (lltl2f 13tl4) t NRO ' 
RIN=R0 „ 



A 
A 
A 
A 



A . 2.40 
.A-^41 
A"^242 
A 243 
244 
245 
246 

249 
250 
251 
252 
253 
254 
255* 
256 
257 
258 ' 
259 
260 
261 
262 
263 
264 
265 
266 
267 
268 
269 
270 
271 
272 
273 
274 
275 



276 
277 
278 
279 
280 
281 
A 282 
A 283 
234 
285 
286 
287 
288 
289 
290 
291 
292 
293 
294 
295 
2,96 
297 
298 
299 



A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 



ERIC 
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143 




GO TO 15 


144 


12 


R I N a R0* ( 1+ A*UM) 


145 




GO TO 15 


146 


13 


RIN=R0»C0SU*DM) 


147 




GO TO 15 


148 


14 


R1N=R0*SIN(A*DM) 


149 


15 


.CONTINUE 


150 




^if (DMH-KTlMc) iVfiStlo 




1 A 
I O 


tc f mo n pt pi an Tn 1<J 


152 




IF (NRO-D 18,17,18, 


153 


17 


RIN'0.0 


154 




GO TO 19 


155 


18' 


RTMH*RTIME/3600. 


156 




RIN=RU«( 1+A*RTMH) 


157 




GO TO 19 


158 


19 


RY(N) 3 RR/Bfc'TA 


v. 


C 
C 


——FORM B-MATRIX 


159 


C 


BE(ll=0.0 


160 




BE(2)=0.0 



161 
162 



163 



164 
^65 
166 
167 
168 
169 
170 
171 



172 
173 
174 



M75 
176 
177 
178 



C 
C 

c 

c 
c 
c 



20 
21 



C 
C 
C 

22 



C 
C 
C 



23 

C 

C 

c 
c 
c 
c 
c 



BE(3)=0.0 

BE(4)±+2.O*D0TMG*TMI/(DX*VM) 

FORM A-HATRU . 

CALL AMTRX • ( AMX ) 



-CHANGE A-DETERMINANT TO POLYNOMIAL FORM 



CALL EIGEN4 ( AMX,AX> 
M=5 
,M1=H 

IF <AX(1>> 22.20^22 
DO 21 JK*i,4 
AX(JK>=AX( JK+i) 
HTRX=3 
.Ml=4 



•COMPUTE THE EIGENVALUE OF THE^A-MATRIX . 



CALL POLRT ( AX»£QF, MTRX ,E IGNV,R00T I , I ER f Mi ) 

MTRX=4 

ML=Mi-i 



-FIND THE LARGEST EIGENVALUE . 



W0=EIGNV(l) 
DO 23 lE*ltML 
WE=E1GNV<IE) 
H0 = AMAX i ( WE t WO ) 



179 
180 



. , CALCULATE THE TIME INCREMENT , SUCH THAT THE TIME INCRE 

ME NT RECIPROCAL TO THE LARGEST tIGEN-VALUE . 

MAKE SURfc THAT THE RANGE IS BETWEEN 0.005 TO 0.05 SECONO , IF 
SMALLER THAN 0.005 SECOND CHANGE THE VALUE TO 0.005 AND IP LAR- 
GER THAN 0.05 SET IT TO 0.05 SECOND . 

0H*A$S< L0/W0)/10.0^jfc * 
DOH*DH 



300 
301 
302 
303 
304 
305 
306 
307 
308 
309 
310 
311 
312 
313 
314 
315 s 
316 
317 
318 
319 
320 - 
321 
322 
323 
324 
325 
326 
327 
328 
329 * 
330 
331 
332 
333 
334 
335 
336 
337 
338 
339 
340 
341 
342 
343 
344 
•A 345 , 
A 346 
347 
348 
349 
350 
351 
A 3?2 
A 353 
3 $4 
355 
356 
357 
358 
359 



A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 

* A 
A 
A 
A 
A 
A 
A 
A 

A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 



X 
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481 



-183 



184 



185 
}86 



187 
188 
189 
190 



191 
192 
193 
194 
195 
196 
197 
198 
199 
200 
201 
202 



203v 
204 

205 
206 
207 
208 
2€9. 
210 
211 
212 
213 



214 



215 



C 

C • 

c 

c 
c 
c 

c 
c 
c 

t a 



c 
c 
I 

24 



C 
C 
C 

c . 

25 



IF (00H.LT.0H1I DH*DH1 . 
• IF (DDH.GUDH2) OH=0H2 

FORM R-><ATRIX . 

CALL RMTRX (RXJ 

~* — FORM H-HATKIX . J 

CALL HftTRX IHHI,- , * - 

\ - MULTIPLY R-MATRIX WITH COLUMfJ MATRIXES AND MULTIPLY 

H-MATRIX WITH COLUMN' VECTOR PH I . 

f CALL GttN ( RX 9 BE ,MTRX,RB) j 
CALL *&XN (HH,PHI,MTRX,HP) ^ j 

{ ..1 FORM NEW PHI-MATRIX BY ADDING THE. TWO COLUHN MATRIX 



26 
27 

C 

28 

C* 

29 

40 

31 

32 

33 

34 

35 

36 

37 



38 



39 



DO 24 IP=1,MTRX 

PHl(iP)*RB<IP)*HPJ,IP) \ 

WRITE (NW,55) N ,DMH ,RY ( N) , PPW , QW ,TF£,TMG t TMOUT, A AF 

IF CDMH-TE) 25,27,27' , , « 

IF MODERATOR TEMPERATURE EXCEEDING 700 DEGREE F ,PRINT A 

WARNING AND GET CUT . 

IF (PHK4J.GE. 1000.0) GO TO 26 
RR-AM*(PHH4)-TM0|*AAF* (PHI ('3)-TF0>*RIN 
R«RR 

DH S DM*DH 
DMH=DH*3600.0 

IF (DMH.GE.TE) GO T€ 27 ^ ... 

N=N«-1 
MX (N)=DMH 
GO TO 7 "* 
4WRITE : (NWt56) 
NN=N 

IF (NOPLT.EO.O) A G0 TD 28 

~ PLOT THE RESULT SIMULTANEOUSLY IN ONE GRAPH . 

CALL TPLOT <MX, MPE, MPEL , HQU , MT EEM, MIEEF/RY, NN ) 
*^TOP 

FORMAT (121 

FORMAT (2F10.5,2F10.7,UO,2F10.5) 
FORMAT (3F10.6) 
FORMAT (3F10.5t2F10.2i2U0) 
FORMAT (2F10.2) 
' FORMAT (3F10.4) > 
FORMAT (3F10. 2) 

FORMAT MxM3H**«********************./tlX,19HM0pULE 2 , FEEDBACK, 
1/, U,23H**i«^^ Ut lO^PUT DATA, /, IX, 14H*** 

^F0RMAT**5X*4HBET A, 5X, 6HLAMBDA, 5X, cui 1 ^/^ *X^i 8 

1L REACTIVITY, /,12X,UH( SEC**-1 1.6*. 7H< SEC >,18X*5H( % l,/,3X,f8 

PERI0D,/.9X,5H<SEC),14 

1X,7H< l/SECI,/,8X,F7.3,13X,F7.3f/) 



360 
361 
362 
363 
364 
365 
366 
367 
368 
369 
370 
371 
372 
A 373 
A 374 
375 
376 
377 
378 
379 
380 
381 
382 
383 
384 
385 
386 
387 
388 
389 » 
390 
391 
392 
393 
394 
395 
396 
397 
398 



A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 
A 

* 
A 

« 

A 



>A 399 



400 
401 
402 
403 
404 
405* 
406 
407 
408 
409 
410 
411 
412 
413 
414 
415 
416 
417 
A M 418 
A*419 



\ 



129 



' - 52 - ? 

216 40 FORMAT < 5X, 14HC0QL ANT COEFF . , 5X, 10HC0NSTANT 8 , 5X, 20HRES0NANCE ESC. 

1 PR0B.,/,6X,12H< , Jt7X,8H< ),7X f 18H< ) 

2 f /,6X,E12.5,7X, E8.2,13X,F5.3,/ ) 

217 41 " FORMAT (5X,11HFUEL RAOI US» 5X, 5HPI TCH,$X, 8HFUEL CP . t 5X, 12^FUEL OENS 

1ITY»5X,11HFUEL HE IGHT ,5X t 1 8H , ,5Xtl2H *t 

2/t7X,6H( FT ),8X,6H< FT )»2X,12H( 8TU/L0-F )»3X,12H< LB/FT**3 ) ,7X 
3t6H( FT >,9X,16H »7X,10H , / , 7X, F7 .5, 7X» F7. 

45t4XfF6.4t8XfF6.2tl0XfF6.2f /I 

218 42 FORMAT t 5X, 18HNUMBER OF ASS EMBLY , 5X , 26HNUMBER OF ROO PER ASSEMBLY, 

l/.5Xt 18HI ) ,5X,26H< )„/.l2X, 

2X4, 23X, I4t/J * ■ - 

219 43 fORMAT ( 5X , 16HC00LANT VELOC ITY, 5X f 20HINLET COOLANT TEMPT .,/, 8X, 10H 

IT FT/SEC ),15X,5H< F I . /. 1 1X.F4. 1 f 18X.F5.1 • /I 
220- 44 FORMAT 1 9X, 6H0PT.ICN, 4X, 4HPL0T , 9X, 1 1 , /, 19X, 1 IHTYPE INSER. » 2X » 1 1 » / ) 

221 45 FORMAT ( 5X» 13HI NI TI AL POWER ,9X,4H< MW) , IX ,E1 1 .4, / ,5X, 19HGUESSE0 RUE 
! * 1L TEMPT. ,4X,3H<F),5X,F6.1,/,5X,22HGUES*E0 COOLANT TEMPT •« 1 X,3Ht F ) , 

26X,F5.1,/// f 2X,17HSN0 OF INPUT OATA . / , 2X, 18 ( 1H*) , ///) 

222 46 FORMAT <//,5XA2H ENO TIME ,3X,14HTIME INCREMENT t /, 8X , F6.3» 10X, F 

16.4,//) 

223 47 FORMAT <5X,22H CONVERGENCE FACTOR * ,F9.5,//> 

224 48 FORMAT (// , 5X , 1 18H*** CHECK DIMENSION t IFIX{TE/0H) HAS TO BE SMA 

1LLER OR EQUAL TO THE DIMENSION OF MPEL , MPb ^MTEtfF , MTEEM, MY,'R Y AND M 
2QU . t///) 

225 49 FORMAT < /// ,5X, 27HEQUI LI BR I UM STATE I NI TI AL LY , / f 4X , 28 ( 1H*) , /// I 

226 50 FORMAT .( 5X,5HPDWER,5X, 9HPRECU$S0R , 5X, 1 1HFUE L TEMPT ., 5X, 14HC00L ANT 

, lTEHPT. t 5X,10HREACTIVITY,5X,riH • . , / f 3X, E10 .4, 2X f E10 .4, 6X, F 

Z7.2, 1 IX ,F6.2, 11 X, F6.3 ,////) , 

227 51 FORMAT ( 85 HH*) »////» 5X » 10H REYNOLDS *, 3X, 9HPR ANOTL #, 3X, 9HC0LBURN 

1#,3X,9HMASS FL0W,3X,1$HHEAT TRANFER COEFF. , / ,42X, 9H< LB/HR ),4X,18 
' 2HTBTU/HR-F-FT**2 I , / ,6X, F 8 . 1 ,6X, F4. 1 , 5X,F 7.3 ,5X,E 1 1. 4,8X ,F7. 1 , /// 
31 * * 

228 52 FORMAT < 5X, 2HND , 10X, 4HTI ME 1 6X ,'10HRE ACTI V ITY ,7X, 5HP0WER 1 1 1X» 9HPRECU 

lRS0R»/t4X«4H( >,8X,5H(S6C),9X,3H($), 11X,4H(MW) ,14X,4H(MW)) 

229 53 FORMAT 17X,9HFUEL TEMP , 5X , 9HM0D. T EMP, 3X,9HEX IT TEMP,8X, 11HFUEL CO 

1EFF.,/,10X,3H(F|,11X,3H(F) ,9X,3H(F) ,14X,5H(1/F) ,//) 

230 54 ' FORMAT < // , 8X , 19H. . . NO FEEDBACK ...,//) 

231 55 FORMAT OX, 1 4 f 8X, F7.3 1 7X ,F5f.2 ,7X, Ell .4, 7X, E 11 .4, / , 8X, F7.2 1 7X, F7 .2t 

15X,F7.2,9X,€11.4) ^ 

232 56, FORMAT <//,5X,39H CRITICAL TEMPERATUR E,//) 

233 / £ND 



420 
4?1 
422 
423 
424 
425 
426. 
427 
428 
429 
430 
431 
432 
433 
434 
435 
436 
437 
438 
439 
440 
441 
442 
443 
444 
445 
446 
447 
448 
449 
450 
451 
452 
453 
454 
455 
456 
457 
458 
459 



234 




SUBROUTINE RMTRX (R) 

• * 


& 


1 


C 


B 


2 




C 




B 


3 




c 


* ^FORM R-MATRI X , DI AGONAL ' MATR I X IN EQUATION 2.4.20 


B 


4 




* c 


SUPPORTING ROUTINE NONE ~ 


B 


5 




c 




B 


6 




C 


\ 


B 


7 


2*5 
236 




DIMENSION R(4t4)«A<4t41 


B 


8 




COMMON BETAfXLtX,FX,RF,DH f FOENStCPF, VF ,FH , VM, OX , OOTMG, WO , HX , RR , C 


B 


9 






1PMX»PI» NROO 


B 


10 


237 




CALL AMTRX <A) 


B 


11 


238 




00 3 1=1,4 


B 


12 


239 




00 3 I I»l f 4 


B 


13 


240 




. IF CI-II) 1#2»I 


B 


14 


241 


1 


R< 1,111 = 0.0 


B 


15 


242 




GO TO 3 . 


B 


• 16 


243 


2 


K(I f 1I) = (EXP( AlI,II)*0H)-1.0l/A<I, 11) 


B 


17 


244 


3 


CONTINUE 


B 


la 



y 



130 



V 



-,53 



245 
246 



247 



248 
249 

250 

251 

252 

253 

254 

255 

256 

257 

258 

259 

260 

261 

262 

263 

264 

2b5 

266 

267 

26*8 



269 



C 

C . 

c 



c 
c 
c 
c 
c 
c 
c 
c 

c , 

c 
c 
c 
c 
c 
c 

. c 
c 
c 
c 
c 
c 
c 
c 
c 
c 



RETURN 
ENO 



SUBROUTINE HMTRX (,H ) % , 

t FORM H- MATRIX AS IN THE EQUATION 2*4.19 

SUPPORTING ROUTINE NONE 



JiS^^t^JuI'.^iF.OH.FOENliCPF.VF.FH.VM.OX.OOTMG.HO.HX.RR.O 

1PMX,PI,NR0D 
CALL AMTRft <A) 
H(Ul) = EXP<AiltU*DH) 

HClt2)**tlt2)*CEXPCMO*OH)-EXPUIltll*DHll/«WO-A(l,lll 
H< lt3) = 0.0 

HI2!l!«AC2til*<EXPCW0*OH)-EXPUI2tZ)*DHII/IW0-Al2t2ll 

Hl2.2) = EXPU<2.2)*DH) 

H<2t3)=0.0 

H!3il)aA(3tl»*IEXPlW0*DH)-EXPU(3t3J*DH))/(WO-A(3',3l.) 
H<3t2)=0.0 

H<3.3MEXP(A(3.3>*DH) - 

H(3,4)=A(3.4)*(EXP<W0*DH)-eXP(A(3,3),*DH)l/(H0-A(3,3)l 
H(4.1^A(4,l)*(EXP(W0*DH)-EXP(A(4,4)*Dti))/<W0-A(4. 6 4)) , 

H<4,2)=0.0 
H<4»3>=0.0 

H<4,4) = EXPUi4,4)*DH> v 

RETURN 

E^D 



SUBROUTINE POLRT ( XCOF ,COF , Mt ROOT R, ROOT 1 1 1 6R, Ml ) 



SUBROUTINE POLRT 



PURPOSE 

• COMPUTES THE REAL AND COMPLEX 



fX Rt) 



10TS Of A REAL POLYNOMIAL 



USAGE - . t ^ 

CALL POLRT (XCGF ,COF,M f ROOTR.ROOTI t IER,Ml) 

DESCRIPTION OF PARAMETERS * _ fAi 

XCGF -VECTOR OF M*l COEFFICIENTS OF THE POLYNOMIAL 

ORDERED FROM SMALLEST TO LARGEST POWER 
COF -WORKING VECTOR Of LENGTH M + l 1 
' M -ORDER-OF POLYNOMIAL ' 1 

ROCTR-RESULT ANT VECTOft^OF LENGTH M CONTAINING REAL RO(tfS 
OF THE POLYNOMIAL > ' , 

" ROOT 1-RESULTANT VECTOR OF LENGTH M CONTAINING THE 

CORRESPONDING *I MAGINARYVROOT.S OF THE POLYNOMIAL 
IER -ERROR C0O£_ WHERE 
1ER=0 1«D ERROR 

M LESS THAN ONE 

M GREATER THAN 36 , 
UNABLE TO DETERMINE ROOJ WITH 500 INJERATIONS 
ON 5 -START ING VALUES 



1ER*1 
IER =2 
IER=3 



;oo IN^E 



ft 

P 


19 


ft 

P 


20 


r 
U 




r 
L 


2 


r 
L 
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r 
L 




r 
L 
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r 
L 




,r 


7 


r 
U 


a 
o 


r 
L 


9 


r 
L 


1 fi 


r 
L 


1 1 

1 1 


L 


1 7 
1 C 


r 
L 




r 
L 




r 
L 




c 


16 


r 


1 7 

L f 


c 


18 


r 


19 


r 
u 


20 


r 


2 1 


r 
L 


7 7 




7 "\ 


£ 

L 


7U 


r 
L 


7** 


r 
L 


c u 


r 
L 


71 
C f 


r 
L 


7 0 


D 


1 
1 


D v 


£ 


V 




V 


A 
*♦ 


r\ 

1/ 


C 
-> 


n 
U 


o 


V 


7 


n 


Q 
O 


V 


Q 




i n 


n 

V 


i i 
1 1 


r\ 
l> 


1 c 


- V 




n 


1 4 


0 


15 


D 


16* 


D 


17 


0 


18 


D 


,19 


0 


20 


D 


21 


D 


22 


, D 


23 


D 


24 


D 


25 


D 


26 



|ER?C 



X31 
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270 
271 



272 



273 
274 
275 
276 
277 



278 
279 



28D 
281 



282 
283 



C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 

c 

c 

G 

c 
c 



1 

c 
c 
c 

2 
3 

c 
c 
c 
4 

c 
c 
c 

5 



Ml 



IER-4 HIGH ORDER COEFFICIENT IS ZERO, 
-NUMBER 0F COEFFICIENT , M+l 

(ADDED ARGUMENT FROM THE ORIGINAL TO GET MORE FLEXIBLE 

'« • DIMENSION ) 



REMARKS 

LIMITED TO 36TH,0KDER POLYNOMIAL DR lESS. 
, FLOATING POINI OVERFLOW MAY OCCUR FOR Hl&H ORDER ' 

POLYNOMIALS 8&T WILL NOT AFFECT THE ACCURACY OF THE RESULTS, 
» ^ J 1 

SUBROUTINES 4ND FUNCTION SUBPROGRAMS REQUIRED 
. NONE , . 

METHOD t * 

NEWTON-*RAPHSDN v ITERATIVE TECHNIQUE. *THE« FINAL ITERATIONS 
DN EACH ROOT ARE PERFORMED USING THE ORIGINAL POLYNOMIAL 
RATHER THAN THE REDUCED POLYNOMIAL TO AVOID ACCUMULATED 
ERRORS IN THE' REDUCED POLVNDMIAC. " 

DIMENSION XCOF(Kl),eOF(Ml) f RDDTRt'M ) , ROOT I < M ) 

DOUBLE PRECISION XC„ Y&, X, Y , XPR |YPR,UX ,UY, V, YT ,XT,U » XT2i YT2i SUMSQ, 
1 DX f OY f TEMP i ALP HA , DAB S 




If . A DOUBLE PRhC.ISIDN VERS I ON OF THIS ROUTINE IS DESIRED, 

C IN COLUMN 1 SHDULO BE REMOVED FRUM THE' DOUBLE PKECISIUN 
STATEMENT WHICH FOLLOWS • 

0 

DOUBLE PRECISION XCCF , CDF, ROOTR, ROOT I 



THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENT 

APPEARING IN OTHSR ROUTINES USED IN CONJUNCTION WITH THIS 
ROUTINE. \ *• 

THE DOUBLE PRECISION VERSION MAY BE MODIFIED BY CHANGING THE 
CONSTANT IN STATEMENT 78 TO 1.0D-12 AND IN STATEMENT 122 TO 
1.00-1D. THJ.S WILL* PROVIDE HIGHER^ PREC I SI DN RESULTS AT THE 
-COST DF EXECUTION TIME 



IFIT=0 

N=M 

IER=0 

IF (XCOF(N+l)) 1,4,1 
IF (N) 2,2,6 



—4 



IERsl 
RETURN 



-SET ERROR CODE* TD 1 



—SET ERROR CODE 'TD 4 



IER-4 
GO TO 3 



SET ERROR CODE TD 2 



IER»2 * 
GO TO 3 
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71 

72 
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1-6 



IF W4-36) 7Wt5.« 

NX*N 

NXX=N*1 

N2*l ' 
KJ1*N*1 

8 L*ltKJl 
MT»KJ1-L*1 
COF<MT) = XCOF(L) 



-SET ItflTIAL VALUES 



X0=. 00500101 
YO-0.01000101 



IN=0 
X«XO 



ZE R0 INITIAL VALUE COUNTER 



-INCREMENT INITIAL VALUES AND COUNTER 



X0=-10.0*YO 
t YO=-10iO*X 



„ X=XO 
Y*YQ, ____ 
I.N«IN+1 
*G0 TO 12 
IFiT«l 
XPR*X 
YPR= Y 



-SET X AND Y TO CURRENT VALUE 



-EVALUATE -POLYNOMIAL AND DERIVATIVES 



ICT»0 

tftaO.O 

UY=0.0 , 

V=0.0 

YT*0.0 

XT=1.0 

U=C0F(N*1> 

IF (U) 14 f 27tl4 

DO 15 I»ltN 

L«N-I*l 

TEMP=COF(LI 

XT2«x*xT-y*rr 

YT2«X*YT*Y*XT 
U=U*TEMP*XT2 
V=V*TEMP*YT2 
F 1 = 1 

UX*UX*FI*XT*TEMlP 

UY*UY-FI*YT.*TEMP 

XT«XT2 

YT*YT2 

SUMSQ*UX*UX*UY*UY 
IF (SUMSO) ffet23tlc> 
OX»<V*UY-U*UXl/SUHSC 
X*X*DX 
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- 56l - 

329 ' DY*-CU*UY*V*UX)/SUMSQ 

330 " Y*Y*DY 

331 IF <OABSCDY)*DABS<DX)-1.0D-5) 21,17,T7 
C 

C -STEP HERAT ION COUNJER * 

332 17 ICTaiet+l 

333 IF ( ICT-5001 13,18, 18 

334 18 1 IF (IFIT) 21,19*21 

335 19 IF ( IN-5) 10,20,20 

' C 

C SET ERROR CODE TO 3 



i 



C 

336 20 IER-i 

337 GO TO 3 

338 21 DO 22 L»1,NXX 

339 NT*KJ1-L*1 

340 TEMP-XCGFCMT) u 

341 XCOF(MT)*COFCLI 

342 22 CDFU) = TEMP 

343 I TEMP*N 

344 N=NX ■ * 

345 NX=ITEMP 

346 IF UFIT) 25,11,25 

347 23 IF i IFIT) 24,10,24 * 
346 . 24 X=XPR . 

349 Y=YPR 

350 25 X IFIT=0 

351 IF <OABS(Y)-i.0D-4*DABSU) ) 28,26,26 

352 26 ALPHA*X«* 
53* SUHSQ?X*X*Y*Y 

354 N=N-2 

355 * GO TO 29 

356 27 X=0.0 * 

357 NX=NX-1 

358 tfXX=NXX-l 

359 28 Y=0.0 
360' SUMSO=0.0 

361 ALPHA=X 

362 N=N- 1 

363 29 C0F(2)=C0F'U)*ALPHA*C0F< W 

364 IF (N.EU.O) GO TO 31 

365 00 30 L=2 , N 

366 30 COF(L*l)=COFtJ.*l)*ALPHA*COF<L)-SUHSQ*CDF(L-l) 

367 31 R.00TICN2I«Y * 

368 R00TR(N2)=X 

369 N2*N2*i 

370 • IF (SUMSO) 32,33»,32 

371 32 Y=-Y * 

372 SUMSQ-0.0 

373 GO TO 31 

374 33 IF (N) 3,3,9 

375 END 



376 SUBROUTINE TPLOT (Ml , M0,M2 . M5,M6, M8, M9, JX) 
C 

C \ 

c TPLOT IS PLOTTING SEVERAL VARIABLES IN ONE 

C DUES NOT REPRESENT ANY VARIABLE , IT IS INTEGER 
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GRAPH . THE X-A 
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SEQUENCES • 


E 


5 



TP 



134 



33 - * 

it,* 
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C 

' c 
c 
c 

377 
376 
379 ' 
360 

382 * 

363 

364 

385 
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39p 
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399 
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40* 
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FOR NEGATIVE' VALUfcS . THE ZERO LINE IS IN THE 0.5 LINE • 
SUPPORTING ROUtlNfc NONE* 

PIHENSION LIN6l6l)tINUMl9) 
INTEGER PL, MI ,S5tBLtSLtS9tS0 f S6tSltS2 
READ <5,8) PL.MI,$5 f BLrSLtS9 f S0,Sl,S^tS6 
MXY=0.0 
MIN0*0.0< 
MIN5-0.0 
MIN6*0.0 
MIN2=0.0 
Ml N8=G.O 
MIN9*0.0 

PHl5*0.0 yf 
PHI6=0.0 
PHI0*0.0 
PHI2=*0.0 
PHI8=0.0 
PHI9=0.0 
00,1 I^ltJX 
IF IMINO.GT.MOI I) ) 
(MlN5.fcT.M5ll)) 
(HIN6.GT.M6l I) ) 
(MIN2.GT.M21 I) ) 
(MIN8.GT.M8II)) 
(MIN9.GT.K9I I )) 
( ABS ( M0( I ) ).GT.PH10) 
m ( ABS ( M2( I ) ).GT.PHI2) 
(ABSIM8U) ).GT,PHI8) 



IF 
. I F 
IF 
IF 
IF 
IF 
IF 
IF 



IF UjBS(M9(I) I.GT.PHI9) 
IF (ABS(M5( I) I.GT.PHI5I 
IF (ABS(M6( I) J.GT.PHI6) 

CONTINUE 
JJ=JX 

J J0*JJ*6*1 
J Jl-JJ* 1 
WRITE (6t9) 
WRITE C6tl0) 
PHlO^PHlO+ABStMlNO) 
PHI5=PHI5*ABS(MIN5> 
PHI6=PHI6*ABS(MIN6) 
PHI2=PHI2*ABS(MIN2) 
PHI8=PHI8*ABS(MIN8) 
PHI9*PHI9*ABS(MIN9) 
00 2 I»1«JJ 
IF (MINO.LT.0.0) 
(MIN5.LT .0.01 
(MIN6.LT.0.0) 
(MIN2.LT.0.0) 
(MINB.LT.O.O) 



MIN0 = M0( I ) 
Ml N5 S M5 ( I ) 
MlN6=M6( I ) 
MIN2=M2( I) 
MlN8*M8( I ) 
MIN9*M9( I) 

PHI 0=ABS ( MO ( I 
PHt2*ABS(H?( I 
PHI8=*BS(M8( I 
PHI9=ABS4M9< I 
PHI5=ABS(M5( I 
PHl6=AftS(M6( I 



IF 
IF 
IF 
IF 

IF (MIN9.LT.0.0) 
M0< I)*MO(I)/PHIO 
M5(I)*M51II/-PHI5 
M6(I)»M6(I|/PH16 
M2U)-M2UI/PHI2 
M8UI=M8(I>/PHI8 
M9(I)«H9(I)/PHI9 



MOU )*MO(I )*ABS<MIN0) 
M5( I ) = M5< I )*ABS(MIN5) 
M6( I)«MfrlI)*ABS«MlN6) 
M2II) = M2U)*AB.$(M!N2) 
M8 ( I )=M8(I)*ABS(MlNfl) 



M9(I)=»M9m+ABS(MlN9) 
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33 
434 * 3 
435 
436, 
437 
436 
439 
440 
441 
442 

443 % 

444 

445 

446 

447 • 

448 

449 

450 

451 * 

452 

453 

454 

455 

456 

457 s 

458 

459 

460 

461 4 

462 

463 

464 

465 

446 

4*7 

468 

469 

470 

471 

472 

473 

474 

475 , 

476 



4*7 

478 

479* 

480 

461 

482 

463 

484 

465 

466 

487 

468 

469 

490 



00 3 I=lt9 
I NUM ( I )~ I 
WRITE (6«1L) 
0(3 7 1=1, JJ1 
IF (I. ECU GO TO 5 
MXY*MU IM) 
IP8«M8( I-l)*60*L0 
IP5=M5( I-l)*60*1.0 
IP9-M91 I-l)*60U.O 
IP6«M6( I-U*60*U0 
IPO*MO( I-i)*60+L0 
IP2=M2( I-1)*60*L0 
00 4 Il = W56,5 
L I NE ( II )»BL 
DO 4 12=1,4 
13=1,1*12 
IF Ul.EQ.IPO) 
IF (I3.EQ.IP0) 
IF UI.EQ.IP2) 
IF (.I3.fO.IP2) 
IF (Il.E0.IP5) 



( INUHl \ \ ,I?l,9) 



LINE(Il)=SO 
L I NE ( t3 ) = SO 
LINE(U)=S2 
L I NE ( I 3) -S2 
LI NE ( I 1 ) = S5 



IF f ( I3.E0.'IP5) t LINt( I 3) =S5 
.1/ (It.E0.IP6) LINEUl) = S6 
IF (I3.EQ.IP6) LINE(I3)=S6 
IF ( I l.EO.IPd) L I NE ( I I ) = SL 
fF (I3.E0.IP8) LINE(I3)=$L 
IF (I1.EQ.IP9) LINE(Il)*S9 
IF CI3.E0.IP9) LINt( I3)=S9 
CONTINUE 
LINE(6l)=PL 
IIL=I-i 

IF ( LEO. IPO) LINE(1)=SG 
IF (I.E0.IP2) _tJNE(l)*S2 
IF ( I.E0.IP5) JLlNFrtJ'SS 
IF ( I .tO. 1^6) L INt ( i')=S6 
IF (1.E0.IP8) L INE( 1 )=SL 
IF '( I.EQ.1P9J LINEII)»S9 
IF (IP0.EQu6l) LINE(6l)=S"C 
(IP5.E0.6l ) 
( IP2.E0.6l) 
( IP6.EQ.6i ) 
( IP8.F.Q.61) 
(IP9.E0.6l) 



LINE(61)=S5 
LINE(6l)=S2 
LlNE(Gl)=S6 
LINE(6l)=SL 
LINE(61)*S9 



IF 
IF 
IF 
IF 
IF 

IF ( IP0.NE.i.0R.£P2.Ni:.i.0R.IP5.NE . I .OR . t P6.NE . I .OR . I P8. NE. I. OR. IP 
19.NE.11 LINE(i)=PL % 
. WRITE (6,12) MXY t (LlNE(KK) , KK = 1 ,6 1 ) 
IF ( I.EO^JJi) GO TCT 7 
CONTINUE 
00 6 Il«l,56,S 

00 6 t2=l,4 . 
13=11+12 . / ' 
LINEU3)=BL 
CONTINUE 

coittTnue , ~ 

' WRITE (6fl3) 

.WRITE (6,15) 

WRITE (6,16) 

WRITE (6,14) 
STOP 



( I NUMt I), 1=1, 9) 

PL tHI tS5 t'BL#SLtS9tS0t SI, S2tS6 
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,491 & 

492 9 

493 *0 
49*, 11 
49ST 12 

4*f /I* 

498 15 

499 16 



500 ' 



FORMAT (HAD 
FORMAT (1H1) 

FORMAT (35X, 16HRELATIVE DENSITY) ■ 

FORMAT U7X,9(2X,2H0*tIltlH )/ 1 14X f 10( 6H* ),1H*) 

FORMAT UX f 4HTlMEttX f F7.3, 1X.61A1) 

FORMAT ( 14X , 10( 6H*~ ) i 1H*. / , 16X i 9< 3 X,2H0 . • 1 1 1 ) 

FORMAT UH1) - ' 

FORMAT (15X.17HINPUT CHARACTER .llAlt/) 

FORMAT (2^H ,/,15X f lOH P -^fOHfeR./ t 15X t 18 

IH L - LOG. OF* POWER f/ ?15X t 22H Q - PRECURSOR DENS ITYrAt 15X ?21H F - 
2FUEL TEMPERATURE t/ t!5X v 26H K - MODERATOR TE MPERATUREt / t 15X , 15H R - 
3 REACTIVITY./, i5Xt42HlF* VARIABLE PLOTTED HAS NEGATIVE VALUE TtfE./t 
415X,41HAXIS FOR THfS VARIABLE I S° SHIFTED TO THE ,/, 15Xt 16HCENTER 0 
5F Y-AXIS,///) 

END 
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SUBROUTINE GXN <G,XN»N,GX) 



MULTIPLY THE G - MATRIX WITH THE INITIAL VECTOR CO 

LUMN TO GET THE NEXT ITERATION . 

SUPPORTING ROUTINE^ NONE 

DIMENSION G(NtN) ? XN ( N ) f GX < ) 

00 2 I»ltN 

GG»0.0/ > 

DO 1 J=lfN 9 < . 

GG=GG+G( It J)*XN( J j 

GXCII=GG 

RETURN 

END ? 



SUBROUTINE AMTRX (A) 



FOR* A-MATRIX IN EQUATION 2.4.12 

SUPPORTING ROUTINE NONE 



I? DIMENSION A(4,4) 

* COMMON BETA.XLfX.FX.RFf DH, FOENS ,CPFt VF.FH, VM, OXt DOTMG, WO.HX, RR , C 
IPMXtPItNROD # ~ ' 

«AU,1>=(RR-BETA)/XL 
A(lt2)=X 
A(l,3)»0.0 

A(lt4)=0.0 * - y * 

,A(2t 1>=BETA/XL t 
At2.2)=-X • 
AC2t3)=0.0 
• A(2t4)=0.0 

A(3»1)=1.0/<2.0*OPF) 
A(3,2)»0.0 1 ♦ 

A(3t3)*-4*C*PI*FX*FH*NR0D/CPF 
A(3.4)= 1 -A(3t3) ' 
A(4 t n = 1.0/(D**<;PHX*VM>, 
A(4t2>*0.0 
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ft ' 
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52 ? 
528 
529 
530 



53 i 
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c 
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c 
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532 
533 
534 
535 



536 

537 
538 
b39 
540 



A<4,3) = 0<,0 

Al4»4i>*-2i0*D0TMG/<DX*VM) 

RETURN 

ENO 

L 

SUBROUTINE EIGEN4 (A,B) 



EIGEN4 IS SUBSTITUTING DETERMINANT A TO POLYNOMIAL FORM 

SUPPORTING ROUTINE NONE 



OIMENSION AI4,4)»B15) 
DOUBLE PRECISION. B 

B(1)=A< L.l)*A12,2)*A<3,3)*A(4*4l-AI 1 , 2 ) *A( 2 1 1 1 *A ( 3,3)*A ( 4, 4 ) 
B(2)=-(A(1, i)*M^rZ ( )*A«3t3) t) +A( i » 1 ) *A ( 2, 2)*A < 4, 4) +A ( i, 1 )*A( 3, 3 )*A ( 4 
l,4)+A«2,2)*A«3,3)*A(4»4)-A«it2)*A(2,l)*A(3»3)-A(i f 2)*A(2,l)*A(4 f 4) 

2) * ^ 

B«3)=A(3»3)*A«4,4) + AU,l)*A(2»2)+AUti)*A(3»3H"A<it i)*A(4 l 4)+A(2»2 

1)*A( 3»3)*-A«2, 2>*A(4,4)-A( i»2)*A(2*l) 
BI4)=-1A( ltil*A12t2)+A(3t3)+AI4,4) ) 
BI5>=1.0 

RETURN P 
END 



G 25 

G 26 

G 27 

G 28 



H 
H 
H 
H 
H 
H 
H 

H 



1 

2 
3 
4 
5 
6 
7 
8 
9 



H 10 
H il 



12 
13 

^ 
15 
16 
•17 
18 
19 



//DATA 



\ 



k 

W" "o 
jERIC 



^0 



138 
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*********************** 

MODULE 2 • FEEDBACK 
*********************** 



INPUT OAT A 
************** 



BETA LAMBDA NEUTRON GEN. TIME INITIAL REACTIVITY 

( SEC**-1 P •* < SEC ) i $ ) 

0,006450 0,07695 v 0.10*)E-03 0,30 

. ■ ■ 

INSERTION TIME CONSTANT PERIOD 

<SEC> U/SEC) 

1,000 1.500 

COOLANT COEFF. CONSTANT B RESONANCE ESC*, PROB, 
C ) C ) * ( ) 
-0.50000E-»04 ' ******** , ^ 0.800 

FUEL RADIUS PITCH FUEL CP. FUEL DENSITY FUEL HEIGHT 

( FT ) i FT ) ( BTU/LB-F ) ( LB/FT**3 ) ( FT ) 

0.01504 0.04733 0.0590 43.20 . 12.00 

NUMBER OF ASSEMBLY - NUMBER <CF ROD PER ASSEMBLY 
( ) ( ) 

145 9 208 

COOLANT VELOCITY, INLET COOLANT TEMPT. 

< FT/SEC ) ( F ) 

13.0- 400.0 

OPTION PLOT 1 
TYPE INSER, V 

INITIAL POWER , (MW) 0.1000E04 

GUESSED FUEL TEMPT. (F) 500.0 « 

GUESSED COOLANT TEMPT. (F) 200.0 



ENO OF INPUT DATA 
****************** 



END TIME TIME INCREMENT 
8.000 0.0100 



CONVERGENCE FACTOR = 0.01000 



. ; ' '• 139. • . 



1ERJC 



EOJILIBRIUJ* STATE INITIALLY 
**************************** 



- PpWER PRECURSOR 
O.fOOOe *04 0.838£E 06 



FUEL TEMPT. 
504.31 



COOLANT TEMPT. 
412.73 



REACTIVITY 
0.000 • 



******************************************* 



REYNOLOS # PRANDTl"* C0L8URN #' MASS FLOW HEAT TRANFER COEFF. 

c ( LB/HR > ( BTU/HR-F-FT**2 ) 

509503.6 1.0 0.04* 0.U72E 09 9044.7 



NO TIME REACTIVITY 

I i , (SEC) > ; (t) 

FUEL T&MP MOD.* TEMP EXIT TEMP 

(Fl (F) (?) 



1 

2 
3 
4 

,5 
6 
7 
8 
9 
10 
11 
12 
13 
14 
15 
1* 



$04.31 
496.42 
497.10 
510.64 
521.65 
529.63 
535.33 
539/41 
542.34 
544.44 
545.97 
547*0^ 
547.91 
548.53 
549.01 



0.000 
0*<f46 
0.096 
0.146 
0.196 
0..246 

0.346 
0.396 
0.446 
0.496 
0.546 
0.596 
0.646 
P. 696 
D.746 



412.73 

413.04' 

412.92 

^13.39 

413.88 

414.32 

414*71 

415.06 

41^.37 

415.65 

415.90 

416.12 

416.31 

4 lti .49 

416.64 



0.00 0, 

• 425. v 47 
0.30 , 0. 

4^6.08 
0.30 ♦ 0, 

425.84 * 
0.30 0, 

426. V* 
0.29 0, 

427.76 
0. ; 29 0 

428.64 
0.29 " , 0 

429.42 
0.28 0 

430.12 
0.28 0 

430.74 
0.28 0 

431.30 
(T.28 0 

451.79 
6.28 0 

' 432. 2i 
0.28 0 

432.62 
0.27 0 

432.97 
0.27 0 

433.29 , 
0.27 0 



POWER 
l^W) 

.1000E 
.9345E 
.1373E 
.1420E 
.14226 
.1419E 
.1416E 
.14J3E 
.1410E 
•1408E 
•1406E 
.1404E 
.14038 
•1403E 
•1402V 
.14'02E 



PRECURSOR 
(MW) 

FUEL COEFF. 
(170) 



04 

0.4968E- 
03 

0.5008E- 
04 

0.5004E- 
04 

0*.4937E- 
04 

0.4885E- 
04 

0.48 t 48E- 
04 

0.4822E- 
04 

0.4804E- 
04 

0.4791E- 
04 

0;4782E- 
04 

0.4775E 
04 

0.4770E 
04 

0.4766E 
04 

0.47^46 
04 

0.4762E 
04 



0.8382E 
06 

0.8381E 
•06 

0.83T9E 
•06 

0.8391E 
•06 

0.8404E 
•06 

0.8418.E 
-06 

0.8431E 
-06 

0.8444E 
-06 

0.845/E 
-06 

0.8470E 
-06 

0.8483E 
-06 

0.8496E 
-06 

0.8508E 
-06 

0.8521E 
-06 

0.8533E 
-06 • 

0.8546E 



06 
06 

0 \ 

06 
06 
06 
06 
06 
06 
06 
06 
06 
06 
06 
-06 
06 
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549*38 




416. 78 




A 2 1 CI 




0.4760E- 


-06 






17 




'0 .796 




0.27 




r\ 1 a rtic 


04 


0.8558E 


uo 






549*69 




4 1 6 0-9 1 








0.4759E- 


-06 






18 




0 . 846 




U • Z f 




n 1 aa?c 
U • 1 *tU£C 


04 


0.8570E 


yj 0 






549>94 




417.02 




A 3 A t\ C 




0.4758E- 


-06 






» 19 


v 


0. 896 




0.27 




f\ 1 a Air 


04 


0.8583E 


06 




550 .17 




417. 13 




434.26 




0.4757E* 


-06 






20 




0.946 




0.27 




0 • 1403E 


04 


0.8595E 


06 






550.37 




417. 22 




434.44 




0.4756E- 


-06 






21 




0 .996 


- 


0.27 




0 • 14036 


04 


0.8607E 


06 






550.56 




41 T .31 




434.61 




0.4755E-06 






22 




1 .046 




0.27 


4, 
f 


0. 1404E 


04 


0.8619B 


06 






550.74 




417.38 




434. 76 




0.4754E 


-06 






23 




1 .096, 




-0.03 




0. 1405E 


04 


0.8631E 


06 






550.91 




417 .45 




434.90 




0.4753E- 


-06 






24 




1. 146 




-0.03 




0 • 101 IE 


04 


0.8643E 


06 






550.95 




417. 52 




435. 03 




0.4753E- 


-06 






25 




1 .196 




-0 .03 




f\ aoonc 
U • WoUc 


03 


0.8643E*06 






539. 15 




417.05 




434. 10 




0.4805E 


-06 






26 




i .246 




-0.03 




a 1 n a n c 
0 . 1 000c 


04 


0.8642E 


06 






530.19 




416.62 




433.23 




0.4845E 


-06 






27 




1 .296 




-0.03 




0. 1 003E 


04 


0.8641E 


06 






5^3.75 




416.23 


* 


432 .46 




0.4875E 


-06 






28 


1. 346 




-0. 02 




0 • 1 UUPfc 


04 


0.8640E 


06 






519.15 




415 .88 




43 1 .77 




0.4897E 


-06 






29 




1 .396 




*0.02 




0.1 007E 


04 


0.8639E 


06 


* 




51 5*87 




415. 58 




43 1 • 1 1> 




0.4912E 


-06 






4 

30 




1 .446 




-0 .02 




0 • 1009c 


04 


». 0.8638E 


06 


**> 




513.54 




415. 30 




43 (7. 60 




0.4923E 


-06 






31 




1 .496 




-0.02 




0. 1 01 IE 


04 


0.8637E 


06 






51 1.87 




415.06 




430.1 1 




0.4931E 


-06 






32 




1 .546 




-0.02 




0. 1013E 


04 


0.8637E 


06 






510.67 




414.84 




429.68 




0.4937E 


-06 






33 




1.596 




-0.01 




0. 1014E 


04 


0.8636E 


06 






509.82 




414.64 




429.29 




0.4941E 


-06 






34 




1.646 




-0.01 




0 • 10 16E 


04 


0.8636E 


06 






509.20 




414.47 




428. 95 




0.4944E 


-06 






35 


1 .696 




-0.01 




0 • 1017E 


04 


0.8635E 


06 






508.76 




414.32 




428.64 




0.4946E 


-06 






36 




L l .746 




-0.01 




0. 1018E 


04 


0.8635E 


06 






508.43 




414.18 




428.37 




• 0.4948E 


-06 






37 




1.796 




-0.01 




0. 1019E 


04 


0.8634E 


06 






508.20 




414- 06 




428. 12 




0.4949E 


-06. 






38 




1.846 




-0.01 




0. 1020E 


04 


0.8634E 


06 






508.02 




413.95 




427.91 




0.4950E 


-06 






£9 




,1.896 




-0.01 




0 • 10216 


04 


0.8634E 


06 




507.89 


413.86 




4£7. 72 




0.4S51E 


-06 






40 




1.946 


-d.01 




U • 1 \JC. 1C 


04 


0.8633E 


06 






507.79 




413*77 




427.55 




0.4951E 


-06 




41 




1.996 




-0.01 




0.1022E 


04 


0.8633E 


06 






507,71 




413.70 




427.39 




0.4952E 


-06 






42 




2.046 




-0.01 




O.1023E 


04 


0.8633E 


06 






507.66v 


413.63 




427.26 




0.4952E 


-06 






43 




2.096 




-0.01 




0.1023E 


04 


0.8632E 


06 






507.61 




413.57 


1 


427.14 




0.49526 


-06 






44 




2.146 




-O.Oi 




0.1023E 


04 


0.8632E 


06 






507.57 




413.52, 




427.03 




0.4952E 


-06 






45 




2.196 




-0.01 




'0.1024E 


04 


0.8632E 


06 






507.54 




413.47 




426.94 




0.4952E 


-06' 






46 




2.?46 


-0.01 




0.10246 


04 


0.8632E 


06 
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3U 1 • 3 A 




413.43 




426*8-5 




0.4953E- 


06 




47 




2«.296 


-0.00 


0.10256 


04 


0.8632E 


06 




S07 SO 




413.39 




426.78 




0.4953E- 


06 


d6 


Aft 




2. 346 


-tf.00 




0. 102 5E 


04 


0.8631E 




*>07 Aft 




413.35 




426. 7 I 




0.4953E-06 




AO. 




3 ^ Q A 

A » J70 


-0.00 




0. 1025E 


04 


0.8651 E 


06 , 




HO'7 4 7 




413.33 




A7A 6*5 




0.4953E- 


06 




3 U 




3 .AAA 


> 


-0.00 




0 • 1025E 


04 


0.8631 E 


06 








413.30 




A 3 a AO 

*T A Om O \J 


0* 1025E 


0.4953E- 


06 








2 .496 




-0.00 


A 3 A ( ¥W 


04 *■ 


' 0.8631E 


06 v 




*i07 A 5 




413.28 




0 1076F" 


0.4953 E- 


0b' * A ' 




3 A 




3 5, AA 
A . 3*tO 


-0.00 




-OA 


0.863]. E 


06 




PU / • *T*T 




413.25 




A3 h S 1 




0 .49^3 E- 


06 








3 RQA 




-0.00 




0 j 1 076F 

DC 




.0.863LE 


06 V 


Of. 






413.24 




A3 A A7 

tfc Oft 1 




0.4953E- 


06 






3 AAA 




-0.00 




0 1076F 


04 




06 




3U I .Ha 




413.22 




A3 a A A 
Ml O . *t*t 




0.4953E- 


C6 




33 




7 AQA 
A . 070 


-0.00 




0.1026E 


04 


0.8631 E 


06 / 




*i »1 7 A 1 




413.20 




426.41 




0.4953 E- 


06 * 




3 O 




3 - 7AA 


-0.00 




0.&026E 


04 


0.8630E 


06 - 




<»< r l7 A 1 




413.19 




426.38 




0.4953E- 


06 




57 




3 7Q A 




-0.00 




0.1026E 


04 


0.8630E 


06 




*>fl7 AO. 
3U 1 • *t U 




413.18 




426.36 




0.4953E- 


06 




CO 

• 3<J 




3 ftA A 




-O'.OC 




0. 1026b 


04 


0.8630E 


06 




*i 07 Art 




413. hi 




426.34 




0.49 53£- 


06 




c.o 




3 ft Q A 




-0.00 




0.1026E 


04 


0.8630E 


06 




cm io 
3U I . J V 




413.16 




426.32 




0.4953E- 


06 








3 QAA 

A. « 7*tO 




-0.00 




0.1026E 


04 


0.8630E 


06 




3U 1 • 37 




413.15 




426.30 




0.4953E- 


06 ' 




ol 




3 QQ A 

A . 770 




-0.00 




0.1026E 


04 


0.8630E 


06 




3U «■* 3 V 




413.14 




426.29 




Q.4953E- 


06 




0£ 




OAA 




-0.00 


Sr0.,1027£ 


04^ 


0.%630E 


06 




30 f • Jo 




413.14 




426.28 




0.4953b- 


■06 




A n 




no a' 

3 . \J 7D 


~ -0.00 




0.1027E 


04 


0.8630E 


06 




K ,"l 7 1Q 
3 J / • D O 




413.13 




426.26 




0.4953E- 


•06 




0** 




5 1 AA 




-0.00 


426.25 


0.1027E 


04 


0.8630E 


06 




30 I m D » 




413.13 






0.4953E- 


•06 




A C 
03 




A 1 OA 
3 • 1 7 O 




-0.00 




0.1027E 


04 * 


0.8629E 


06 




'<*fl7 ft 
3U UJO 




413.12 




426.24. 


0*4953E- 


•06 




oo 




3 AA 
3 . A*r O 


-0.00 




0.1027E 


04 


. 0.8629E 


06 




cm *a7 




413.12 


-o.ob 


426.24 




0.4953E- 


-06 




£7 

O f 




^ 3QA 
3 . A 7 O 






0.1027E 


04 


0.8629E 


06 




CA7 57 




413.11 




426.23; 




0.4953E-06 




A Q 

OO 








-0.00 




0.102 7* 


04 


.0.8629E 


06 




507.37 




413.11 




426.22 




0.4953E- 


-06 




« AO 




^ ^ Q A 
3 . J7 O 




-0.0Q 




■Q ? 1027E 


04 


'0.8629E 


06 - 




Cn 7 3. 7 




413.11 




426.2 2 




0.4953£- 


-06 




70 




3.446 




-0.00 




0.10~27E 


04 A 


0.8629E 


06 




507.36 




413.11 




426.21 




0.4953E^06 . 




71 




3-496 




-0.00 




0.1027E 04 


0.8629E 


06 




507.36 




413.10 




426.21 




0.4953E- 


-06 




72 




3.546 




-0.00 




0.1027E 


04 


0.8629E 


06 




507*36 




419.10 




426.20 




0.4.953E- 


-06 1 




73 




3*596 




-0 .00 




0.1027E 


04 


0.8629E 


06 




507.36 




413.10 




426.20 




0.4953E-06 




74 




3.646 




-0.00 




0\1027E 


04 


0.8629E 


'06 




507.36 




413.10 




426.20 




0.4953E-06 




75 


3.696 




-0.00 




0.1027E 


04 


0.8629E 


;06 




507.36 




413.10 




426. IP 




0.4953E-06 




76 




3.746 
% 




-0.00^ 


0.1027E«04 


0.8628E 


06 ' 
9 






1 507.35 


77 






507.35 


78 






507.35 


79 






507.35 


80 






507.35 


81 






507.35 


82 






507.35 


83 






507.34 


84 






507.34 


85 




, 507*34 


86 






507.34 


87 






507.34 


88 






507.34 


89 






507.34 


90 






507.34 


91 






'507.34 


92 






507^33 


93 






507.33 


94 






507.33 


95 




, 507*33 


96 




507.33 


97 






507.33 


98 






507.33 


99 






507.33 


100 






507.33 


101 






507.32 


102 






507.32 


103 






507.32 


10.4 






507.32 


105 






507.32 


106, 



413.10 426.19 0,49536-06 

3.796 -0.00 ' 0.1027E 04 1 0.8628E 06 

413.09 426.19 0.49536-06 

3.846 -0.00 * 0.1027E 04 0.8628E 06 

• 413.09 426.19 0.49536-06 

3.896 -0.00 0.1027?^ 0.8628E 06 

413.09 426.U <JV£9536-06 

3.946 -0.00 0.1027E 04 S> 0.8628E 06 

413.09 ' 426.18 0.49536-06 

3.996 ^0.00 0.1027E 04 0.8628E 06 

41:3.09 426.18 0.49536-06 

4.046 . -0.00 0.1027E 04 0.8628E 06 

413.09 426.18 0.4953E-06 

4.096 -0.00 0.1027E 04 0.8628E 06 

413.09 426. 18 0. 49536-06 

4.146 -0.00 , 0.10276 -04 0.86286 06 

413.09 426.18 0.49536-06 

4.196 <-0.00 , 0.10276 04 0.86286 06 

413.09 - 426.17 * 0.49536-06 

4.246 -0.00 * 0.1027& 04 0.86286 06 

413.0$ 426.17 0.49536-06 

4.296 -0.00 0.10276 04 0.86276 06 

413.09 426.17 4 0,49536-06 . 

4.346 -0.00 4 0.10276 <*4 0.8627E 06 

413.09 426.17 0.49536-06 

4.396 -0.00' 0.10276 04 0.86276.06 

413.09 426.17 • 0.49536-06 t 

4.446 -0.00 0.10276 04 / 0.86276 06 

' • 413.09 426.17 0.49536-06 

4.496 -0.00 0.10276 04 , 0.86276 06 

413.09 ' 426.17 ^ 0.49536-06 

4.546 -0.00 0.10276 04 \ 0.86276.06 

413.09 426.17 0.49536-06 

4,596 -d.00 0.10276 04 , 0.86276 06 

413.08 426. 17 0.49536-06 

4.646 -0.00 ,0.10276 04 0.86276 06 

> 413.08 426.17 0.49536-06 

4.696 -0.00 0.10276 04 0.86276 06 

413U08 " 426. 17 0.49536-06 

4.746 -0.00 0.1027.6 Q4 0.86276 06 

413.08 426.17 0.49536-06 

4.796 -0.00 0.10276 04 0.86276 06 

413.08 , 426.17 0.49536-06 
4.846 -0.00 0.1027E 04 0.86276 06 

413.08 426.17 0.49536-06 

4.896 -0.00 0.10276 04 0.86266 fib 

413.08 426. 17 0.49536-06 

4.946 -0;00 0.10276 Q4 - <t>. 86266. 06 

413.08 426.17 0.49536-06 

4.996 -0.00. 0.1027E 04 0.86266 06 

413.08 426.17 0.49536-06 

5.046 -0.00 0.10276 04 0.86266 06 

413.08 426.17 0.49536-06 

5.096" -0.00 0.10276 04 0.86266 06 

413.08 426.17 0.49536-06 4 

5.1*46* ' -0.00 * 0.10276 04 0.86266 06 

413.08 426.17 * 0.49536-06 

5.196 -Q.00 0.10276 04 0.86266 06 

> 413.08 426.17' 0. 49536-06 ' 

5.246 -0.00 0.10276 04 0.86266 06 
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